A parameter uniform almost first order convergent numerical method for a non-linear system of singularly perturbed differential equations

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 differential 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 normally 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 derived. 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 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 interaction 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 catalysts.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 where k 1 , k −1 and k 2 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 with the initial conditions s(0) = s 0 , e(0) = e 0 , 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 Often the enzyme concentration e 0 is very small compared with that of substrate concentration s 0 ; in such situations it follows that e 0 s 0 = ε 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. where It is assumed that the nonlinear terms satisfy for all f defined on [0, 1] × R 2 .These conditions ( 5), ( 6) and the implicit function theorem [8], ensure that a unique solution ), exists for problem (4).
For any vector-valued function y on [0, 1] the following norms are introduced: i=0 is a real valued function defined on Ω N .The discrete maximum norm for such mesh functions is defined by 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 The reduced problem corresponding to ( 7) is given by Since r ∈ C 2 , f (x, r 1 , r 2 ) 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, and the singular component w(x) satistfies The bounds of the derivatives of smooth component are given in As in [3] and [5] the smooth component v is further decomposed as follows: v = q + q.
Using (14) in (15) then gives Consider now the linear operator This operator satisfies a maximum principle, see [1].

III. SHISHKIN MESH
For the case ε 1 < ε 2 , a piecewise uniform Shishkin mesh Ω N with N mesh-intervals is constructed on Ω = [0, 1] as follows.The interval where Then, on the subinterval (τ 2 , 1] a uniform mesh with N  2 meshpoints is placed and on each of the sub-intervals (0, τ 1 ] and (τ 1 , τ 2 ], a uniform mesh of N 4 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 using the backward Euler scheme on the piecewise uniform mesh Ω N .The discrete problem is where 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 Proof: where and ξ ik (x j ), η ik (x j ), 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, 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 N 0 > 0 such that, for all N ≥ N 0 , U − u ≤ CN −1 lnN.
Here, C and N 0 are independent of ε 1 , ε 2 and N. Proof: Let x ∈ Ω N .From the above lemma, 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 (C 6 H 12 O 6 ) n is the substrate, amylase is the enzyme and maltose n(C 12 H 22 O 11 ) is the product, which is represented as follows: 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 Here x, u 1 , u 2 are the dimensionless quantities as in [6] and are given by, x = k 1 e 0 t, u 1 (x) = c(t) e0 , u 2 = 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 u 1 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 u 2 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 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].The notations D N , p N and C N p * denote the ε-uniform maximum pointwise two-mesh differences, the ε-uniform order of convergence and the ε-uniform error constant respectively and given by D N = max 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, it could be noted that as the domain increases u 1 and u 2 approach zero since, as the time increases the concentrations s and c become zero.In both the graphs, the component u 1 , exhibits an initial layer whereas u 2 exhibits no layer.In Table 1, it is seen that the parameter uniform order of convergence p N is monotonically increasing towards the value 1 as N increases.This is in agreement with Theorem 1.
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 financial support through the Rajiv Gandhi National Fellowship to carryout this work.

D N N p * 1 − 2
−p * .Then the parameter -uniform order of convergence and error constant are given by p * = min N p N and C * p * = max N C N p * .

TABLE 1
Values of D N ε , D N , p N , p * and C N p * for α = 0.9.