www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE A parameter uniform almost first order convergent numerical method for a non-linear system of singularly perturbed differential equations R.Ishwariya∗, J.Princy Merlin†, J.J.H.Miller‡, S.Valarmathi§ ∗Department of Mathematics, Bishop Heber College, Tiruchirappalli, TamilNadu, India. ishrosey@gmail.com †Department of Chemistry, Bishop Heber College, Tiruchirappalli, TamilNadu, India. pmej 68@yahoo.co.in ‡Institute of Numerical Computation and Analysis, Dublin, Ireland. jm@incaireland.org §Department of Mathematics, Bishop Heber College, Tiruchirappalli, TamilNadu, India. valarmathi07@gmail.com Received: 30 September 2015, accepted: 11 August 2016, published: 11 September 2016 Abstract—In this paper, a biochemical reaction namely Michaelis-Menten kinetics has been modeled as an Initial Value Problem (IVP) for a system of singularly perturbed first order nonlinear differen- tial equations with prescribed initial values. A new numerical method has been suggested to solve the problem of Michaelis-Menten kinetics. The novelty is that, a variant of the continuation technique suggested in [2] with the classical finite difference scheme is used on an appropriate Shishkin mesh instead of the usual adaptive mesh approach nor- mally used for such problems. In addition to that, a new theoretical result is included which establishes the parameter uniform convergence of the method. From the computational results inferences are de- rived. Keywords-Singular perturbation problems, boundary layers, nonlinear differential equations, finite difference schemes, Shishkin mesh, parameter uniform convergence. I. INTRODUCTION A differential equation in which small parameters multiply the highest order derivative and some or none of the lower order derivatives is known as a singularly perturbed differential equation. In this paper, an initial value problem for a sytem of singularly perturbed nonlinear differential equations is considered. In [9], a parameter uniform numerical method for a class of singularly perturbed nonlinear scalar initial value problems is constructed. Also results for a problem where two reduced solutions intersect are discussed. In [11], the asymptotic behaviour of the solution for a nonlinear singular singularly perturbed initial value problem is studied. In [3], a numerical method, for a system of singularly perturbed semilinear reaction-diffusion equations, involving an appropriate layer-adapted piecewise Citation: R.Ishwariya, J.Princy Merlin, J.J.H.Miller, S.Valarmathi, A parameter uniform almost first order convergent numerical method for a non-linear system of singularly perturbed differential equations, Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 1 of 8 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... uniform mesh, is constructed and it is proved to be parameter uniform convergent. In [10], a numerical method for a singular singularly perturbed initial value problem for a system of nonlinear differential equations is suggested. Singularly perturbed differential equations have wide applications in Biology. In [5], a parameter-uniform numerical method for a model of tumour growth, which is a first order semilinear system of singularly perturbed delay differential equations, is suggested and is proved to be almost first order convergent. One of the applications of singularly perturbed differential equations in the field of Biochemistry is the modelling of enzyme kinetics. In [12], analytical approximations to the equation for modelling Michaelis-Menten kinetics are provided. In [13], the enzyme kinetic reaction scheme originally proposed by V. Henri is considered and a method of computing the rate constants, whenever time course experimental data are available is also demonstrated. In [14], a new approach is developed for the study of the Michaelis-Menten kinetic equations at all times t, based on their solution in the limit of large t. In the present paper, an enzyme kinetic reaction is considered in the form of an initial value problem for a system of singularly perturbed nonlinear differential equations of first order. The changes in the solution components are studied and illustrated numerically. It is also proved that the numerical method considered is parameter uniform and is almost first order convergent. A biochemical reaction is a process of interac- tion of two or more substances to produce another substance. Biochemical reactions are continually taking place in all living organisms and in all body processes. Such reactions involve proteins called enzymes, which act as remarkably efficient cata- lysts. Enzyme kinetics is the study of the chemical reactions that are catalysed by enzymes. In enzyme kinetics, the reaction rate is measured and the effects of varying the conditions of the reaction are investigated. Enzymes react selectively on definite compounds called substrates. One of the most basic enzymatic reactions, that occurs in Henri- Michaelis-Menten kinetics (1900) and Michaelis- Menten kinetics (1913), involves a substrate S reacting with an enzyme E to form a complex SE which in turn is converted into a product P and the enzyme. We represent this schematically by S + E k1 k−1 SE k2→ P + E, (1) where k1,k−1 and k2 are constant parameters associated with the rates of reaction. As in [6], by using the Law of Mass Action in (1), a system of nonlinear differential equations is obtained ds dt = −k1es + k−1c, de dt = −k1es + (k−1 + k2)c, (2) dc dt = k1es− (k−1 + k2)c, dp dt = k2c, with the initial conditions s(0) = s0, e(0) = e0, c(0) = 0, p(0) = 0, where s,e,c,p denote the concentration of S, E, SE, P respectively. It is not hard to get p(t) from the last expression of (2). Observing that de dt + dc dt = 0 and introducing the dimensionless quantities as in [6], (2) reduces to du dτ = −u + (u + K −λ)v, ε dv dτ = u− (u + K)v, (3) u(0) = 1, v(0) = 0. Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 2 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... Often the enzyme concentration e0 is very small compared with that of substrate concentration s0; in such situations it follows that e0 s0 = ε � 1. This system is said to be partially singularly perturbed, since only the second equation of the system is singularly perturbed. It is to be noted that the right hand sides of (3) are nonlinear. Motivated by this, the following more general nonlinear system of singularly perturbed differential equations is now considered. ~T~u(x) = E~u ′(x) + ~f(x,u1,u2) = ~0 on (0, 1], ~u(0) = ~u0, (4) where x ∈ [0, 1], ~u(x) = (u1(x),u2(x))T and ~f(x,u1,u2) = (f1(x,u1,u2),f2(x,u1,u2)) T . E is the 2 × 2 matrix, E = diag(~ε ), ~ε = (ε1,ε2) with 0 < ε1 ≤ ε2 ≤ 1. It is assumed that the nonlinear terms satisfy ∂fk ∂uk ≥ β > 0, ∂fk ∂uj ≤ 0, k,j = 1, 2, k 6= j, (5) min 1≤i≤2   2∑ j=1 ∂fi ∂uj   ≥ α > 0 (6) for all ~f defined on [0, 1] ×R2. These conditions (5), (6) and the implicit function theorem [8], ensure that a unique solution ~u ∈ C2 = C × C, where C = C0([0, 1]) ∩ C2((0, 1]), exists for problem (4). For any vector-valued function ~y on [0, 1] the following norms are introduced: ‖ ~y ‖= max{‖ y1 ‖, ‖ y2 ‖} and ‖ yi ‖= sup x∈[0,1] |yi(x)|. The mesh Ω N = {xi}Ni=0 is the set of points satisfying 0 = x0 < x1 < ... < xN = 1. A mesh function V = {V (xi)}Ni=0 is a real valued function defined on Ω N . The discrete maximum norm for such mesh functions is defined by ‖ V ‖ Ω N = maxi=0,1,...,N|V (xi)| and ‖ ~V ‖ Ω N = max{‖ V1 ‖ΩN ,‖ V2 ‖ ΩN} where the vector-valued mesh function is defined by ~V = (V1,V2) T = {V1(xi),V2(xi)}T , i = 0, 1, ...,N. 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 (4) can be rewritten in the form ε1u ′ 1(x) + f1(x,u1,u2) = 0, ε2u ′ 2(x) + f2(x,u1,u2) = 0, x ∈ (0, 1], ~u(0) = ~u0. (7) The reduced problem corresponding to (7) is given by f1(x,r1,r2) = 0, f2(x,r1,r2) = 0, x ∈ (0, 1]. (8) Since ~r ∈ C2, ~f(x,r1,r2) is sufficiently smooth. Also from (8) and the conditions (5) and (6), the implicit function theorem [8] ensures the existence of a unique solution for (8). 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, 1]. A 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) + ~f(x,v1,v2) = ~0, x ∈ (0, 1], ~v(0) = ~r(0), and the singular component ~w(x) satistfies E~w ′(x) + ~f(x,v1 + w1,v2 + w2) −~f(x,v1,v2) = ~0, x ∈ (0, 1], ~w(0) = ~u(0) −~v(0). (9) The bounds of the derivatives of smooth compo- nent are given in Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 3 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... Lemma 1: The smooth component ~v(x) satis- fies |v(k)i (x)| ≤ C, i = 1, 2; k = 0, 1 and |v′′i (x)| ≤ Cε −1 i , i = 1, 2. Proof: As in [3] and [5] the smooth component ~v is further decomposed as follows: ~v = ~̃q + ~̂q. Here ~̂q is the solution of f1(x, q̂1, q̂2) = 0, (10) ε2 dq̂2 dx + f2(x, q̂1, q̂2) = 0, x ∈ (0, 1], (11) with q̂2(0) = v2(0), and ~̃q is the solution of ε1 dq̃1 dx + f1(x, q̃1 + q̂1, q̃2 + q̂2) −f1(x, q̂1, q̂2) = −ε1 dq̂1 dx , (12) ε2 dq̃2 dx + f2(x, q̃1 + q̂1, q̃2 + q̂2) −f2(x, q̂1, q̂2) = 0, (13) x ∈ (0, 1], with q̃1(0) = q̃2(0) = 0. Using (8), (10) and (11), a11(x)(q̂1 −r1) + a12(x)(q̂2 −r2) = 0, (14) ε2 d dx (q̂2 −r2) + a21(x)(q̂1 −r1) + a22(x)(q̂2 −r2) = −ε2 dr2 dx , (15) where aij(x) = ∂fi ∂uj (x,ξi(x),ηi(x)), i,j = 1, 2, and ξi(x), ηi(x) are intermediate values. Using (14) in (15) then gives ε2 d dx (q̂2 −r2) + ( a22(x) − a12(x)a21(x) a11(x) ) (q̂2 −r2) = −ε2 dr2 dx . Consider now the linear operator Lz(x) := ε2z ′(x) + ( a22(x) − a12(x)a21(x) a11(x) ) z(x) = −ε2 dr2 dx , where z = q̂2 −r2. This operator satisfies a maximum principle, see [1]. Thus, ‖ q̂2 −r2 ‖≤ Cε2 and ‖ d(q̂2 −r2) dx ‖≤ C. Using this in (14), ‖ q̂1 −r1 ‖≤ Cε2. Differentiating (14), we get ‖ d(q̂1 −r1) dx ‖≤ C. Hence, ‖ q̂2 ‖≤ C, ‖ dq̂2 dx ‖≤ C and ‖ q̂1 ‖≤ C. Differentiating (15), ε2 d2 dx2 (q̂2 −r2) + a′21(x)(q̂1 −r1) + a21(x) d dx (q̂1 −r1) + a′22(x)(q̂2 −r2) + a22(x) d dx (q̂2 −r2) = −ε2 d2r2 dx2 . Hence, ‖ d2q̂2 dx2 ‖≤ Cε−12 . Differentiating (14) twice and using the above estimates, it follows that ‖ d2q̂1 dx2 ‖≤ Cε−12 . By using the assumption that ε1 ≤ ε2, we find that ‖ d2q̂1 dx2 ‖≤ Cε−11 . Expressions (12), (13) and the mean-value theorem for ~f, lead to ε1 dq̃1 dx + a∗11(x)q̃1 + a ∗ 12(x)q̃2 = −ε1 dq̂1 dx , (16) ε2 dq̃2 dx + a∗21(x)q̃1 + a ∗ 22(x)q̃2 = 0, (17) where q̃1(0) = q̃2(0) = 0. (18) Here, a∗ij(x) = ∂fi ∂uj (x,ζi(x),χi(x)), i,j = 1, 2; ζi(x), χi(x) are intermediate values. Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 4 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... From equations (16)-(18), by using maximum principle [4], we obtain ‖ q̃i ‖≤ C, i = 1, 2, ‖ dq̃i dx ‖≤ C, i = 1, 2, ‖ d2q̃i dx2 ‖≤ Cε−1i , i = 1, 2. Hence the bounds for ~v can be obtained using the bounds of ~̃q and ~̂q. Lemma 2: The singular component ~w(x) satisfies, for any x ∈ [0, 1], |wi(x)| ≤ CB2(x), i = 1, 2, |w′1(x)| ≤ C(ε −1 1 B1(x) + ε −1 2 B2(x)), |w′2(x)| ≤ Cε −1 2 B2(x), |w′′i (x)| ≤ Cε −1 i (ε −1 1 B1(x) + ε −1 2 B2(x)), i = 1, 2, where Bi(x) = e (−αx εi ) . Proof: From equations (9), ε1w ′ 1(x) + s11(x)w1(x) + s12(x)w2(x) = 0, (19) ε2w ′ 2(x) + s21(x)w1(x) + s22(x)w2(x) = 0, (20) for x ∈ (0, 1], and w1(0) = u1(0) −v1(0), w2(0) = u2(0) −v2(0). Here, sij(x) = ∂fi ∂uj (x,λi(x),θi(x)), i = 1, 2; λi(x), θi(x) are intermediate values. From (19) and (20), the bounds of the singular component ~w can be derived as in [4]. III. SHISHKIN MESH For the case ε1 < ε2, a piecewise uniform Shishkin mesh Ω N with N mesh-intervals is con- structed on Ω = [0, 1] as follows. The interval [0, 1] is subdivided into 3 sub-intervals [0,τ1] ∪ (τ1,τ2] ∪ (τ2, 1], where τ2 = min { 1 2 , ε2 α lnN } and τ1 = min {τ2 2 , ε1 α lnN } . Clearly 0 < τ1 < τ2 ≤ 12. Then, on the subinterval (τ2, 1] a uniform mesh with N2 meshpoints is placed and on each of the sub-intervals (0,τ1] and (τ1,τ2], a uniform mesh of N4 points is placed. 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, 1]. In the case, ε1 = ε2 a simple construction with just one parameter τ is sufficient. IV. DISCRETE PROBLEM The initial value problem (7) is discretised us- ing the backward Euler scheme on the piecewise uniform mesh Ω N . The discrete problem is ~TN ~U(xj) : = ED −~U(xj) + ~f(xj,U1(xj),U2(xj)) = ~0, j = 1, ...,N, ~U(0) = ~u(0), (21) where D−~U(xj) = ~U(xj) − ~U(xj−1) xj −xj−1 . This discrete problem (21) is to be solved on the Shishkin mesh defined above, using the continuation algorithm [2]. Lemma 3: For any mesh functions ~Y and ~Z with ~Y (0) = ~Z(0), ‖ ~Y − ~Z ‖≤ C ‖ ~TN ~Y − ~TN ~Z ‖ . Proof: ~TN ~Y (xj) − ~TN ~Z(xj) = ED−~Y (xj) + ~f(xj,Y1(xj),Y2(xj)) −ED−~Z(xj) − ~f(xj,Z1(xj),Z2(xj)) = ED−(~Y − ~Z)(xj) + J(~f, ~u)(~Y − ~Z)(xj) = (~T ′N )( ~Y − ~Z)(xj), Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 5 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... where ~T ′N is the Frechet derivative of ~TN, J(~f, ~u) =   ∂f1 ∂u1 (xj,ξ11,η11) ∂f1 ∂u2 (xj,ξ12,η12) ∂f2 ∂u1 (xj,ξ21,η21) ∂f2 ∂u2 (xj,ξ22,η22)   and ξik(xj),ηik(xj), i,k = 1, 2 are points occuring in the mean value theorem. Since ~T ′N is linear, it satisfies the discrete maximum principle and the discrete stability result [4]. 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 problem (4) and ~U be the solution of the discrete problem (21). Then, there exists an N0 > 0 such that, for all N ≥ N0, ‖ ~U −~u ‖≤ CN−1lnN. Here, C and N0 are independent of ε1,ε2 and N. Proof: Let x ∈ ΩN. From the above lemma, ‖ (~U −~u)(xj) ‖≤ C ‖ (~TN ~U − ~TN~u)(xj) ‖. Since ‖ ~TN~u (xj) ‖=‖ (~TN~u− ~TN ~U)(xj) ‖, it follows that ‖ (~TN~u− ~TN ~U)(xj) ‖ = ‖ ~TN~u(xj) ‖ = ‖ (~TN~u− ~T~u)(xj) ‖ = E ‖ (D−~u−~u ′)(xj) ‖ ≤ E ‖ (D−~v −~v ′)(xj) ‖ + E ‖ (D−~w − ~w ′)(xj) ‖ . Since the bounds for ~v and ~w are the same as in [4], the required result follows. V. NUMERICAL RESULTS The numerical method proposed above is illustrated through the example presented in this section. Experimental data for the hydrolysis of starch by amylase is used in the numerical illustration. In this experiment starch (C6H12O6)n is the substrate, amylase is the enzyme and maltose n(C12H22O11) is the product, which is represented as follows: (C6H12O6)n Starch + nH2O Amylase −→ n(C12H22O11) Maltose By the procedure explained in the introduction, the enzyme kinetics of the above reaction leads to the formulation of the following Example: Consider the initial value problem ε du1 dx = u2 − (u2 + K)u1, du2 dx = −u2 + (u2 + K −λ)u1, for x ∈ [0, 1], u1(0) = 0, u2(0) = 1. Here x,u1,u2 are the dimensionless quantities as in [6] and are given by, x = k1e0t, u1(x) = c(t) e0 , u2 = s(t) s0 . Based on the experiments and analysis carried out in the departments of Zoology and Chemistry of Bishop Heber College, Tiruchirappalli, Tamilnadu, India, and [7], the values of K and λ are taken to be 0.0477mg−1 and 0.0454mg−1. Observations: The solution profile of the above problem exhibits interesting facts. The solution component u1 of ~u, representing the ratio of concentration of the complex at time t to initial concentration of the substrate, exhibits initial layer in the neighbourhood of t = 0. The second component u2 exhibits no layer in the domain of the definition of the problem. The ratio of the concentration of the complex to that of enzyme has a rapid increase initially but Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 6 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... as time increases, this ratio remains smooth and reaches equilibrium position. On the other hand, the ratio of the concentration of the substrate at time t to the initial concentration of the substrate decreases smoothly. The maximum pointwise errors and the rate of convergence for this IVP are calculated using a two mesh algorithm for a vector problem which is a variant of the one found in [2] for a scalar problem. The results are presented in Table 1 for the domain [0,1]. TABLE 1 Values of DN~ε ,D N,pN,p∗ and CNp∗ for α = 0.9. ε Number of mesh points N 512 1024 2048 4096 2−1 0.361E-03 0.181E-03 0.905E-04 0.453E-04 2−3 0.143E-02 0.717E-03 0.359E-03 0.180E-03 2−5 0.217E-02 0.123E-02 0.688E-03 0.380E-03 2−7 0.217E-02 0.123E-02 0.688E-03 0.379E-03 2−9 0.217E-02 0.123E-02 0.688E-03 0.379E-03 DN 0.217E-02 0.123E-02 0.688E-03 0.380E-03 pN 0.819E+00 0.840E+00 0.859E+00 CNp∗ 0.833E+00 0.833E+00 0.821E+00 0.799E+00 Computed order of ~ε−uniform convergence, p∗ = 0.8194415 Computed ~ε−uniform error constant, C∗p∗ = 0.8329233 The notations DN, pN and CNp∗ denote the ~ε-uniform maximum pointwise two-mesh differences, the ~ε-uniform order of convergence and the ~ε-uniform error constant respectively and given by DN = max ~ε DN~ε where DN~ε = ‖ ~U N ~ε − ~U 2N ~ε ‖ΩN , p N = log2 DN D2N and CNp∗ = DNNp ∗ 1 − 2−p∗ . Then the parameter - uniform order of convergence and error constant are given by p∗ = min N pN and C∗p∗ = max N CNp∗. For ε = 2−5 and N = 512, the solution of the example is presented in the form of a figure in the domain [0,1](Figure 1). And from Figure 2, Figure 1. Numerical solution for ε = 2−5 and N = 512 in the domain [0,1]. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 u1 u2 Figure 2. Numerical solution for ε = 2−5 and N = 512 in the domain [0,30]. 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 u1 u2 it could be noted that as the domain increases u1 and u2 approach zero since, as the time increases the concentrations s and c become zero. In both the graphs, the component u1, exhibits an initial layer whereas u2 exhibits no layer. In Table 1, it is seen that the parameter uniform order of convergence pN is monotonically increasing towards the value 1 as N increases. This is in agreement with Theorem 1. ACKNOWLEDGMENT The authors are grateful to Prof.Svetoslav Markov for suggesting improvements to this paper. Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 7 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 R.Ishwariya et al., A parameter uniform almost first order convergent numerical method ... Further, the authors thank the unknown referees for their suggestions and comments which have improved the paper to its present form. The first author sincerely thanks the University Grants Commission, New Delhi, India, for the fi- nancial support through the Rajiv Gandhi National Fellowship to carryout this work. REFERENCES [1] J.J.H.Miller, E.O’Riordan, G.I.Shishkin, Fitted Numeri- cal 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 & Mathematical Computation (Eds. R.J.Knops & K.W.Morton), Champman & Hall/CRC Press (2000). [3] J.L.Gracia, F.J.Lisbona, M.Madaune-Tort, E.O’Riordan, A Coupled system of singularly perturbed semilinear reaction-diffusion equations, A.F. Hegarty (Ed.), et al., BAIL 2008-Boundary and Interior Layers, Lecture Notes in Computational Science and Engineer, vol. 69 (2009), pp. 163172 [4] S.Valarmathi, J.J.H.Miller, A parameter uniform finite difference method for a singularly perturbed linear dy- namical systems, International Journal of Numerical Analysis and Modelling, 7(3), 535-548 (2010). [5] N.Shivaranjani, J.J.H.Miller, S.Valarmathi, A Parame- ter Uniform Almost First Order Convergent Numerical Method for a Semi-Linear System of Singularly Per- turbed Delay Differential Equations, Biomath 3 (2014), 1411041. [6] J.D.Murray, Mathematical Biology: An Introduction (Third Edition), Springer (2002.) [7] Nikola Sakac̆, Milan Sak-Bosnar, Potentiometric Study of α−Amylase kinetics using a Platinum Redox Sensor, International Journal of Electrochemical Science, 7(2012) 3008 - 3017. [8] Gerald B.Folland, Advanced Calculus,Prentice Hall, 2002. [9] Eugene ORiordan, Jason Quinn, Parameter-uniform nu- merical methods for some singularly perturbed nonlinear initial value problems, Springer (2012) 61:579611. [10] R. E. OMalley, JR.andJ. E.Flaherty, Analytical and Numerical methods for Singular Singularly-Perturbed Ini- tial Value Problems, Society for Industrial and Applied Mathematics, Vol. 38, No. 2, April 1980. [11] Mo Jiaqi and Lin Wantao, Asymptotic Solution of Singularly Perturbed Problem for a Nonlinear Singular Equation, Appl.Math.J.Chinese Univ.Ser.B, 2004,19(2) : 187-190. [12] Marko Golicnik, Exact and Approximate Solutions for the Decades-old Michaelis-Menten Equation: Progress- curve Analysis Through Integrated Rate Equations, Bio- chemistry and Molecular Biology Education, Vol. 39, No. 2, pp. 117125, 2011. [13] Stanko Dimitrov, Gergana Velikova, Venko Beschkov, Svetoslav Markov, On the Numerical Computation of Enzyme Kinetic Parameters, Biomath Communications 1 (2014). [14] M. Kosmas, E. M. Papamichael, E. O. Bakalis, An Iterative Solution of the MichaelisMenten Equations, Match Commun. Math. Comput. Chem. 70 (2013) 971- 986,ISSN 0340 - 6253. Biomath 5 (2016), 1608111, http://dx.doi.org/10.11145/j.biomath.2016.08.111 Page 8 of 8 http://dx.doi.org/10.11145/j.biomath.2016.08.111 Introduction ANALYTICAL RESULTS SHISHKIN MESH DISCRETE PROBLEM NUMERICAL RESULTS References