www.biomathforum.org/biomath/index.php/biomath

ORIGINAL ARTICLE

Robust numerical method for a singularly
perturbed problem arising in the modelling of

enzyme kinetics
John J. H. Miller, Eugene O’Riordan

1School of Mathematics, Trinity College
Dublin 2, Ireland

jmiller@tcd.ie
2School of Mathematical Sciences, Dublin City University

Dublin 9, Ireland
eugene.oriordan@dcu.ie

Received: 19 June 2020, accepted: 22 August 2020, published: 12 September 2020

Abstract— A system of two coupled nonlinear
initial value equations, arising in the mathematical
modelling of enzyme kinetics, is examined. The
system is singularly perturbed and one of the
components will contain steep gradients. A priori
parameter explicit bounds on the two components
are established. A numerical method incorporating
a specially constructed piecewise-uniform mesh is
used to generate numerical approximations, which
are shown to converge pointwise to the continuous
solution irrespective of the size of the singular
perturbation parameter. Numerical results are pre-
sented to illustrate the computational performance
of the numerical method. The numerical method is
also remarkably simple to implement.

Keywords-Enzyme-substrate dynamics, nonlinear
system, Shishkin mesh, parameter-uniform conver-
gence.

I. INTRODUCTION

The Henri-Michaelis-Menten system of non-
linear differential equations arises in the mathe-
matical modelling of enzyme-substrate dynamics,
see, for example, [1], [3], [5], [7]. As analyt-
ical solutions are not available, it is necessary
to solve this system numerically. This may be
difficult when the system is singularly perturbed
[1]. Asymptotic expansions associated with this
singularly perturbed problem are discussed in [9],
[7], [8]. Here, to establish a parameter-uniform
pointwise error bound on the numerical approx-
imations, we construct a Shishkin decomposition
[6] for the solutions. This can be viewed as an
alternative to an asymptotic expansion. Based on
this decomposition, an efficient finite difference
method is constructed, which uses a specially

Copyright: c©2020 Miller et al. This article is distributed under the terms of the Creative Commons Attribution License
(CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and
source are credited.
Citation: John J. H. Miller, Eugene O’Riordan, Robust numerical method for a singularly perturbed
problem arising in the modelling of enzyme kinetics, Biomath 9 (2020), 2008227,
http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 1 of 12

http://www.biomathforum.org/biomath/index.php/biomath
https://creativecommons.org/licenses/by/4.0/
https://creativecommons.org/licenses/by/4.0/
http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

constructed piecewise uniform mesh and an ap-
propriate standard finite difference operator to deal
with the steep gradients that appear initially, see
for example [6]. This numerical method is shown
to be parameter-uniformly convergent, in the sense
that its numerical solutions converge to the exact
solution, uniformly with respect to the singular
perturbation parameter, with essentially first order.
A Shishkin mesh was also used in [4] to solve
the same system, but it should be noted that the
hypotheses [4, (5) and (6)], required by the theory
in [4], are not fulfilled by the system considered
in the present paper.

Asymptotic expansions yield useful approxima-
tions when the singular perturbation parameter ε
is sufficiently small for terms of a certain order
O(εn) to be negligible. However, the accuracy
in parameter-uniform approximations is valid ir-
respective of the size of ε and depends only
on the number of mesh points N used in the
computations. Moreover, there can be a debate
[7], [8] about how to choose the dimensionless
variables when using the quasi-steady-state as-
sumption [8], which becomes irrelevant if one
has a parameter-uniform numerical method. The
Shishkin-mesh method developed below does not
require an asymptotic expansion to compute an ap-
proximation and, hence, the method is simpler than
other approaches which generate individual terms
(or approximations) in an asymptotic expansion
(e.g., [9]).

The structure of the paper is as follows. In the
next section the continuous problem is formulated.
In §3 the problem is discretised and the error anal-
ysis is performed in §4. The main results are stated
in Theorems 5 and 6. The paper concludes with
§5 in which some numerical experiments illustrate
the form of the solution, and its initial layer, and
also support the theoretical error analysis.

II. CONTINUOUS PROBLEM

In the basic model of enzyme reactions (e.g.
[7]), a substrate S reacts with an enzyme E to form
a complex SE which is converted into a product
P and the enzyme. The concentrations of these

variables vary with time τ and are denoted here
by lower case letters

s(τ) = [S], e(τ) = [E], c(τ) = [SE], p(τ) = [P ].

These reactions can be modelled by the follow-
ing system of four first order equations for four
unknowns with given initial conditions

ds

dτ
= −k1es + k−1c, (2.1a)

de

dτ
= −k1es + (k−1 + k2)c, (2.1b)

dc

dτ
= k1es− (k−1 + k2)c, (2.1c)

dp

dτ
= k2c, (2.1d)

s(0) =s0, e(0) =e0, c(0) = 0, p(0) = 0; (2.1e)

where the parameters k−1,k1,k2 are reaction rate
constants. Note that

p(τ) = k2

∫ τ
s=0

c(s) ds and e(τ) = −c(τ) + e0.

Hence, we have only two unknown variables (s
and c) to determine. As in [7], we introduce the
following scalings and nondimensionless variables
u and v:

u(t) :=
s(t)

s(0)
, v(t) :=

c(t)

e(0)
, t := (k1e(0))τ

α :=
k2

k1s(0)
, K =

k−1 + k2
k1s(0)

,ε :=
e(0)

s(0)
.

This leads to the following autonomous system of
two coupled nonlinear initial value equations ([7,
eq. (6.13)]):

Find (u(t),v(t)) ∈ C∞(0,T), 0 ≤ t ≤ T such
that, for all 0 < ε ≤ 1,

u′ = −u + (u + K−α)v, t > 0; (2.2a)
εv′ = u− (u + K)v, t > 0; (2.2b)
u(0) = 1; v(0) = 0;K > α > 0. (2.2c)

We observe the following facts

u′(0) = −1, εv′(0) = 1, (2.3)
εu′′(0) = 1 + K−α + ε,
εv′(t) + αv(t) = −u′(t);

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 2 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

and so there will be an initial layer in v and a
weak initial layer in u (as εu′′(0) = O(1) and
u′(0) = O(1)). The exact solution of problem (2.2)
is unknown.

Lemma 1. For all t > 0, the solutions of (2.2) are
bounded as follows:

0 < u(t) < 1 and 0 < v(t) ≤
1

1 + K
. (2.4)

Proof: Given the initial conditions (2.3) and
the fact that u′,v′ are continuous, there exists an
interval (0, t1), where u′(t) < 0 and v′(t) > 0 for
all t ∈ (0, t1). Hence, there also must exist some
0 < t0 < t1 such that

0 < u(t) < 1 and

0 < v(t) <
u(t)

u(t) + K
< 1, t ∈ (0, t0).

Also, for any positive time t > t1, where u(t) +
K > 0,

v′(t) > 0, if v(t) <
u(t)

u(t) + K
and

v′(t) < 0, if v(t) >
u(t)

u(t) + K
.

Assume that there exists a t∗ > t1 such that

0 < v(t∗) =
u(t∗)

u(t∗) + K
< u(t∗) < 1,

where v′(t∗) = 0 and u′(t∗) < 0. If no such t∗ is
reached, then we are done. If this time t∗ exists,
then v will have a maximum at this point and so
v(t) < 1 for all time where u(t) + K > 0 .

Hence we have that u(t∗) and v(t∗) are both
positive. By (2.2a), there does not exist a least time
t2 > t∗ where u(t2) = 0,u′(t2) ≤ 0 and v(t2) >
0. Moreover, by (2.2b), there does not exist a least
time t2 > t∗ where v(t2) = 0,v′(t2) ≤ 0 and
u(t2) > 0. Hence, if either u or v were to become
negative, then u,v would have to simultaneously
become zero at t2. Hence, by (2.2a) and (2.2b),
we would have that u(t2) = u′(t2) = v(t2) =
v′(t2) = 0. Moreover, by differentiating (2.2a) and
(2.2b), we would find that all derivatives of u (and
v) were zero at this point t2. For a smooth function
u ∈ C∞(0,T], it cannot be that all derivatives

of a non-trivial function at a certain time t2 are
zero. Hence, this point t2 does not exist. It follows
that u(t) > 0, v(t) > 0, ∀t ≥ 0. Finally, as
u(t) + K−α > 0,

u′ = −u + (u + K)(1 −
α

u + K
)v

≤ −u + u(1 −
α

u + K
)

= −
uα

u + K
< 0.

Hence u(t) ≤ 1. Note that g′(z) > 0 for g(z) =
z

z+K. Hence, v(t) ≤ v(t∗) ≤
1

1+K.
From the bounds in this Lemma, we can deduce

that, for i = 1, 2,∥∥∥diu
dti

∥∥∥ ≤ C(1 + ε1−i), ∥∥∥div
dti

∥∥∥ ≤ Cε−i;
where ‖g‖ := maxt∈[0,T ] |g(t)|. However, these
bounds are not sufficiently sharp for our purposes,
as they do not identify the fact that the large
derivatives will only occur initially. To generate
sharper bounds, we will construct a Shishkin de-
composition [6] of the solution.

From ([7, eq. (6.26)]) and by formally setting
ε = 0 in (2.2) and ignoring v(0) = 0, the reduced
solution (u0,v0) is given by

v0(t) =
u0(t)

u0(t) + K
, u0(t) +K ln u0(t) = 1−αt.

(2.5)
The reduced solution (u0,v0) or outer solution
approximates the solution (u,v) outside a neigh-
bourhood of t = 0. Using the stretched variable
τ = t/ε, the solution (u,v) is approximated ([7,
eq. (6.31)]) initially by the inner solution (uI,vI)

uI(t) = 1, vI(t) =
1

1 + K
(1 −e−(1+K)t/ε),

for t ≤ Cε. This motivates the following Shishkin
decomposition of the solution (u,v).

Lemma 2. The solutions of problem (2.2) can be
decomposed as follows:

u(t) = u0(t) + Ru(t), (2.6a)

where
du0
dt

= −
αu0

u0 + K
, u0(0) = 1;

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 3 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

v(t) =
u(t)

u(t) + K
−B(t) + Rv(t), (2.6b)

B(t) :=
1

1 + K
e−

∫
t

s=0

u(s)+K
ε

ds;

and the remainder terms Ru,Rv are bounded by

‖Ru‖ + ‖Rv‖≤ Cε, (2.6c)∥∥∥diRu
dti

∥∥∥ + ∥∥∥diRv
dti

∥∥∥ ≤ C(1 + ε1−i), i = 1, 2.(2.6d)
Proof: The function v(t) can be written in the

form

v(t) =
u(t)

u(t) + K
−

1

1 + K
e−

∫
t

s=0

u(s)+K
ε

ds + Rv(t).

By inserting this expansion for v into (2.2b), we
see that the remainder term Rv satisfies the initial
value problem

ε
dRv
dt

+(u+K)Rv = ε
d

dt

( u(t)
u(t)+K

)
(2.7)

=
εK

(u+K)2
du

dt
; Rv(0) = 0.

Hence,

Rv(t) =

∫ t
s=0

g(t)e−
∫
t

r=s

u(r)+K
ε

drds,

where

g(t) =
K((u + K−α)v −u)

(u + K)2
.

Note that ‖g‖ ≤ C and by the previous lemma,
u(t) + K≥K > 0. Then we deduce that∣∣Rv(t)∣∣ ≤ C ∫ t

s=0
e−
K(t−s)

ε ds ≤ Cε.

Using the differential equation in (2.7), we con-
clude that∥∥∥diRv

dti

∥∥∥ ≤ C(1 + ε1−i), i = 1, 2.
We also have the following decomposition

u(t) = u0(t) + Ru(t),

where
du0
dt

= −
αu0

u0 + K
, u0(0) = 1.

Note that u0(t) is implicitly defined in (2.5). By
inserting the expansions for u and v into (2.2a),
we see that

u′ = −u+(u+K−α)v
= −u+u−

αu

u+K
+(u+K−α)(Rv−B)

= u′0 +
αu0
u0 +K

−
αu

u+K
+(u+K−α)(Rv−B),

where

R′u =
αu0

u0 + K
−

αu

u + K
+ (u + K−α)(Rv −B)

=
αK(u0 −u)

(u0 + K)(u + K)
+ (u + K−α)(Rv −B).

Hence, the remainder term Ru satisfies the initial
value problem

dRu
dt

+
αK

(u + K)(u0 + K)
Ru

= (u + K−α)(Rv −
1

1 + K
e−

∫
t

s=0

u(s)+K
ε

ds);

Ru(0) = 0.

Hence, as (u0 +K)(u +K) ≥K2 > 0, by writing
out a closed form representation for the function
Ru, we deduce that

‖Ru‖≤ Cε and
∥∥∥diRu
dti

∥∥∥ ≤ C(1 + ε1−i), i = 1, 2.
III. DISCRETE PROBLEM

Consider the following implicit linear finite dif-
ference scheme for the continuous problem (2.2):
For all tj ∈ ΩN , find (U(tj),V (tj)) such that: For
all tj > 0

D−U(tj) + (1 −V (tj−1))U(tj)
−(K−α)V (tj) = 0, (3.8a)

εD−V (tj)−U(tj)+(U(tj−1)+K)V (tj) = 0, (3.8b)

and
U(0) = 1; V (0) = 0; (3.8c)

where

D−U(tj) :=
(
U(tj)−U(tj−1)

)
/kj, kj :=tj−tj−1

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 4 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

and the piecewise uniform mesh ΩN is constructed
below. Note that for any s > 0, since |u0(t)| ≤ 1,

u0(s) = 1 +

∫ s
0

du0
dt

dt

= 1 −α
∫ s

0

u0
u0 + K

dt > 1 −
α

1 + K
s.

Hence, the initial layer function B(t) is bounded
by

e−
∫
σ

s=0

u(s)+K
ε

ds ≤ Ce−
∫
σ

s=0

u0(s)+K
ε

ds

≤ Ce−
σ

ε
(1+K− ασ

2(1+K)
)
.

Based on the decomposition (2.6b) and (2.6a), we
will use a piecewise uniform Shishkin mesh [6],
denoted here by

ΩN := {tj|0 ≤ tj := tj−1 + kj ≤ T},

where there the fine and coarse mesh sizes k,K
are defined by

kj = k :=
2σ

N
, j ≤

N

2
;

kj = K :=
2(T −σ)

N
, j >

N

2
.

The transition point σ is taken to be

σ := min{0.5T,
ε

0.5 + K
ln N}.

For simplicity of exposition we assume 1that
ε(ln N)2 ≤ C then σ2 ≤ Cε and for N sufficiently
large,

e−
∫
σ

s=0

u(s)+K
ε

ds ≤ Ce−
∫
σ

s=0

u0(s)+K
ε

ds

≤ Ce
ασ2

2ε(1+K) e−
σ(1+K)

ε ≤ CN−1.

Remark 3. If ε ≤ CN−1 then for t ≥ σ we have

v(t) =
u0(t)

u0(t) + K
+ CN−1

and
u(t) = u0(t) + CN

−1.

In other words, outside the initial layer, the so-
lutions (u,v) are computationally close to the

1If ε(ln N)2 > C then a classical argument can be used to
deduce the error bound in Corollary 7.

reduced solutions (u0,v0) when ε is sufficiently
small.

We next establish that, within the fine mesh, the
discrete solutions U,V are bounded by the same
bounds as their continuous counterparts u,v. In the
next section, we will extend this result to the mesh
points outside the fine mesh. These bounds on
the discrete solutions ensure that the linear system
(3.8) has a unique solution.

Lemma 4. For the solution of (3.8) and N suffi-
ciently large (independent of ε), we have, for all
0 < tj ≤ σ,

0 < U(tj) < 1 and 0 < V (tj) <
1

1 + K
.

Proof: By explicitly solving the linear system
(3.8) at the first internal mesh point, we see that

0 < V (t1) =
k/ε

1 + k + k 1+K
ε

+ k
2

ε
(1 + α)

<
k

ε

1

1 + K
<

1

1 + K

0 < U(t1) =
1 + k(1 + K)/ε

1 + k + k 1+K
ε

+ k
2

ε
(1 + α)

< 1.

Define the associated system matrix

Mj :=I+kj

( 1−V (tj−1) −(K−α)
−1
ε

U(tj−1)+K
ε

)
; (3.9)

and write the discrete problem (3.8) in the form

Mj

(
U(tj))
V (tj)

)
=
(
U(tj−1))
V (tj−1)

)
,(

U(t0)
V (t0)

)
=
(

1
0

)
.

Note that

det(Mj) = 1+kj(1−V (tj−1))+kj
U(tj−1)+K

ε

+
k2j
ε

(
U(tj−1)+α−V (tj−1)(U(tj−1)+K)

)
= 1+kj(1−V (tj−1))+

Kkj
ε

(1−kj)+
kj
ε
U(tj−1)

+
k2j
ε

(
(1−V (tj−1))(U(tj−1)+K)+α

)
Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 5 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

and
M−1j =

A

det(Mj)
, (3.10)

where

A:=I+kj

( U(tj−1)+K
ε

K−α
1
ε

1 −V (tj−1)

)
.

We now complete the argument by using in-
duction. Assume the statement is true for all
1 ≤ i ≤ j − 1. If 0 < V (tj−1) < 1 and
U(tj−1) > 0 then det(Mj) > 1 and M

−1
j > 0.

Hence, U(tj) > 0,V (tj) > 0. Moreover, we can
rewrite

det(Mj) = 1+kj(1−V (tj−1))+
kj
ε
U(tj−1)(1+K)

+
k2j
ε

(
(1 + α−V (tj−1)(U(tj−1) + K)

)
+

(K−kj)kj
ε

(1 −U(tj−1)).

If (1 + K)V (tj−1) < 1 and U(tj−1) < 1 then
V (tj−1)(U(tj−1) + K) < 1 and

det(Mj) > 1+kj(1−V (tj−1))+
kj
ε
U(tj−1)(1+K).

From this, we deduce that

V (tj) =
1

det(Mj)
(
(1+kj(1−V (tj−1)))V (tj−1)

+
kj
ε
U(tj−1)

)
<

1

(1 + K)det(Mj)
(
(1 + kj(1 −V (tj−1)))

+
kj
ε
U(tj−1)(1 + K)

)
<

1

1 + K
.

We rewrite

det(Mj) = 1+kj
U(tj−1)+K

ε
+kj(K−α)V (tj−1)

+kj
(
1−(1+K)V (tj−1)+V (tj−1)(α−

kj
ε

(K−α))

+
kj
ε

(U(tj−1) + α)(1 −V (tj−1)
))
.

If (1 +K)V (tj−1) < 1 and j ≤ N/2 then kj ≤
CεN−1 ln N and, under these assumptions,

det(Mj)>1+kj
(U(tj−1)+K)

ε
+kj(K−α)V (tj−1).

Hence

U(tj) =
1

det(Mj)
(
1 +

kj
ε

(U(tj−1) + K)

+kj(K−α)V (tj−1)
)
< 1.

IV. ERROR ANALYSIS

Consider the error (U − u,V − v)(tj) at each
mesh point tj, which satisfies

D−(U−u)+(1−V (tj−1))(U−u)−(K−α)(V −v)
= u′ −D−u + (V (tj−1) −v(tj))u;

εD−(V −v)−(U−u)+(U(tj−1)+K)(V −v)
= ε(v′ −D−v) + (u(tj) −U(tj−1))v.

We rewrite these equations in matrix form

Mj

(
(U −u)(tj)
(V −v)(tj)

)
=

( kj(u′ −D−u + (V (tj−1) −v(tj))u)
kj
ε

(
ε(v′ −D−v) + (u(tj) −U(tj−1))v

) ),
where the matrix Mj is defined in (3.9). The next
theorem establishes an error estimate within the
fine mesh region [0,σ].

Theorem 5. For the solutions of (2.2), (3.8) and
N sufficiently large (independent of ε), we have,
for all tj ≤ σ,

|(U −u)(tj)| + |(V −v)(tj)| ≤ CN−1(ln N)2.

Proof: For tj ≤ σ,

|
(
u′ −D−u

)
(tj)| ≤ Ck‖u′′‖≤ CN−1 ln N,

|u(tj) −u(tj−1)| ≤ k‖u′‖≤ CN−1;
ε|(v′ −D−v)(tj)| ≤ Cεk‖v′′‖≤ CN−1 ln N,
|v(tj) −v(tj−1)‖ ≤ k‖v′‖≤ CN−1 ln N.

Hence,

|u′ −D−u + (V (tj−1) −v(tj))u| ≤
CN−1 ln N + C|(V −v)(tj−1)|;

|ε(v′ −D−v) + (u(tj) −U(tj−1))v| ≤
CN−1 ln N + C|(U −u)(tj−1)|.

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 6 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

As in the previous lemma, det(Mj) > 1 and using
the explicit form of the inverse M−1j in (3.10), we
deduce that

|(U −u)(tj)| ≤ Ck
(
|(V −v)(tj−1)

+CN−1 ln N(1 + |(U −u)(tj−1)|)
)
;

|(V −v)(tj)| ≤ C
k

ε

(
|(U −u)(tj−1)

+CN−1 ln N(1 + |(V −v)(tj−1)|)
)
.

Hence, using an iterative argument, for all j ≤
N/2

|(U−u)(tj)|+|(V −v)(tj)| ≤ Cj(
k

ε
)2

≤ CN(N−1 ln N)2

≤ CN−1(ln N)2.

It is not as straightforward to establish this error
bound in the coarse mesh region. Before we ex-
amine this error, we present a discrete analogue to
the solution decomposition established in Lemma
2. We first construct the decomposition in the fine
mesh region, as it is relatively straightforward to
do so. Then, in the next section, we inductively
generate this same discrete decomposition and
simultaneously establish error bounds at each point
tj in the coarse mesh region.

We have the following decomposition of the
discrete solution V :

V (tj) =
U(tj)

U(tj)+K
−

1

1+K
BN (tj)+R

N
v (tj),

(4.11a)
where

εD−BN (tj)+(U(tj−1)+K)BN (tj) = 0,
BN (0) = 1. (4.11b)

The remainder term RNv satisfies R
N
v (0) = 0 and

the finite difference equation

εD−RNv + (U(tj−1) + K)R
N
v

=−εD−
U(tj)

U(tj)+K
+kj

U(tj)D
−U(tj)

U(tj)+K
(4.11c)

=
(−εK(U(tj−1)+K)−1 +kjU(tj))

U(tj)+K
D−U(tj).

For tj ≤ σ and N sufficiently large,

U(tj−1) + K = u0(tj−1) + K±CN−1

> 1 −
α

1 + K
σ + K±CN−1

≥ 0.5 + K.

Then in the fine mesh(U(tj−1) + K
ε

)
k ≥

2 ln N

N

and from [6, Lemma 5.1]

Π
N/2
j=1

(
1 +

(U(tj−1) + K
ε

)
k
)−1 ≤ CN−1.

Hence, we have the bound

0 < BN (σ) ≤ CN−1. (4.11d)

We next establish a bound on (RNv −Rv)(tj) for
all tj:

εD−(Rv−RNv )(tj)+(U(tj−1)+K)(Rv−R
N
v )(tj)

= εD−Rv(tj) + (U(tj−1) + K)Rv(tj)+

εKD−U(tj)
(U(tj−1) + K)(U(tj) + K)

−
kjU(tj))

U(tj) + K
D−U(tj)

=ε(D−Rv−R′v)(tj)+(U(tj−1)−u(tj))Rv(tj)+G(tj),

where

G(tj) :=

εK
(

D−U(tj)

(U(tj−1)+K)(U(tj)+K)
−

u′(tj)

(u(tj)+K)2

)
−

kjU(tj)

U(tj) + K
D−U(tj).

Also

D−U(tj)

(U(tj−1) + K)(U(tj) + K)
−

u′(tj)

(u(tj) + K)2

= T(u(tj) + K)2D−U(tj)
−Tu′(tj)(U(tj−1) + K)(U(tj) + K)

= T(u(tj) + K)2(D−U(tj) −u′(tj)) +
Tu′(tj)((u−U)(tj))((u(tj) + K) +
Tu′(tj)(U(tj) + K)(u(tj) −U(tj−1),

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 7 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

where

T :=
1

(U(tj−1) + K)(U(tj) + K)(u(tj) + K)2
.

Note that

u′(tj) = (K−α)v(tj)−(1−v(tj))u(tj)
D−U(tj) = (K−α)V (tj)−(1−V (tj−1))U(tj)
|u(tj)−u(tj−1)|≤Ckj‖u′‖.

If

0 < U(tk),V (tk) < C, for tk = tj−1, tj

then

|G(tj)| ≤
C(ε+ N−1) max

k=j−1,j
{|(U−u)(tk)|, |(V −v)(tk)|}

+CN−1. (4.12)

Hence, using a discrete comparison principle 2

|(RNv −Rv)(tj)| ≤
C(N−1+max

k<j
|(U−u)(tk)|+max

k≤j
|G(tk)|). (4.13)

Since U(tj) > 0, by Lemma 3, we have that for
all tj ≤ σ

|(U −u)(tj)| + |(V −v)(tj)| ≤ CN−1(ln N)2

and hence

|(RNv −Rv)(tj)| ≤ CN
−1(ln N)2, tj≤σ. (4.14)

Theorem 6. For the solution of (3.8) and N
sufficiently large (independent of ε), we have, for
all tj ≥ σ,

0 < U(tj) < 1, 0 < V (tj) <
1

1 + K
and the following parameter-uniform error bound

|(U −u)(tj)| + |(V −v)(tj)| ≤ CN−1(ln N)2.

Proof: The proof is by induction. By the pre-
vious two lemmas the statement is true for tk = σ.

2 If a(tj ) > 0 and W (tj ) is any mesh function such that
εD−W (tj ) + a(tj )W (tj ) ≥ 0,∀tj > 0 and W (0) ≥ 0 then
W (tj ) ≥ 0,∀tj ≥ 0.

Assume now that it is true for all σ ≤ tk ≤ tj−1.
The argument in Lemma 2 to establish the bounds

0 < U(tj), 0 < V (tj) <
1

1 + K
is still valid within the coarse mesh. However, the
argument to establish the upper bound U(tj) < 1
requires an alternative argument in the coarse mesh
region.

Under the induction assumption, for all σ ≤
tk ≤ tj the decomposition in (4.11a)

V (tk) =
U(tk)

U(tk) + K
−

1

1 + K
BN (tk) + R

N
v (tk),

where

εD−BN (tk) =−(U(tk−1)+K)BN (tk), BN (0) = 1

is still applicable. Moreover, as U(tk−1) +K > 0,
we see that

0 < BN (tk) < B
N (σ) ≤ CN−1

and

|(RNv −Rv)(tk−1)| ≤ CN
−1, for all k ≤ j.

Returning to the end of the proof of Lemma 2, we
have that

det(Mj) = 1 + kj

(
U(tj−1) + K

ε

)
+kj(K−α)V (tj−1)+kj

(
1−(1+K−α)V (tj−1)

)
+
k2j
ε

(
U(tj−1) + α− (U(tj−1) + K)V (tj−1)

)
.

Using the decomposition (4.11a), we see that

U(tj−1) + α− (U(tj−1) + K)V (tj−1)
> α±CN−1 > 0

for N sufficiently large. Hence, as in the proof of
Lemma 2, U(tj) < 1.

For each tj, using the decompositions (2.6b) and
(4.11a) we can write the error in v in the form

(v −V )(tj) =
u(tj)

u(tj) + K
−

U(tj)

U(tj) + K
+CN−1 + (Rv −RNv )(tj)

=
K(u(tj) −U(tj))

(u(tj) + K)(U(tj) + K)
+CN−1 + (Rv −RNv )(tj). (4.15)

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 8 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

For tj > σ, using (2.6b) we have that

ε(v′−D−v)(tj) =
1

kj

∫ kj
s=kj−1

ε(v′(tj)−v′(s))ds

=
1

kj

∫ kj
s=kj−1

(u−(u+K)v)(tj)−(u−(u+K)v(s))ds

=
(u + K)(tj)

kj

∫ kj
s=kj−1

v(tj) −v(s) ds + CN−1.

Hence, using (2.6b) again, we get that

ε|(v′ −D−v)| ≤ CN−1.

Also

|u(tj) −u(tj−1)| ≤ CN−1|u′| ≤ CN−1

and for tj > σ

|v(tj) −v(tj−1)| ≤∣∣ u(tj)
u(tj) + K

−
u(tj−1)

u(tj−1) + K
∣∣ + CN−1

≤ CN−1.

In the coarse mesh

|
(
u′−D−u

)
(tj)| ≤ CN−1 + |

(
R′u−D

−Ru
)
(tj)|.

The remainder term Ru can be further decomposed
into the sum

Ru = Su + w,

where
‖S′′u‖≤ C,‖w

′‖≤ Ce−
Kt
ε .

Hence, on the coarse mesh,

|
(
u′ −D−u

)
(tj)| ≤ CN−1.

Hence, the error (U −u)(tj) satisfies

D−(U −u) + (1 −V (tj−1))(U −u)

−
K(K−α)

(u(tj) + K)(U(tj) + K)
(U −u)

= (V −v)(tj−1)u + (K−α)(RNv −Rv)(tj)

+CN−1 + C
kj
ε
e−
Ktj−1
ε .

By the induction hypothesis and using the expres-
sions in (4.13) and (4.12), we deduce that

D−(U−u)+
(
1−v(tj)−

K(K−α)
(u(tj)+K)(U(tj)+K)

± C(ε + N−1)
)
(U −u)(tj)

= (V (tj−1)−v(tj))(U−u)(tj)+CN−1(ln N)2;

D−(U −u) +
( K(U(tj) + α)

(u(tj) + K)(U(tj) + K)
± C(ε + N−1)

)
(U −u)(tj)

= CN−1(ln N)2.

Since the coefficient of the zero order term in
this error equation is strictly bounded below by
a positive number, we deduce that

|(U −u)(tj)| ≤ CN−1(ln N)2.

The bound on the error in V (tj) follows from
(4.15). The proof is completed by induction.

Let Ū and V̄ denote the piecewise linear inter-
polants of U and V , respectively. Then, as in [6],
we can extend the result to a global error bound.

Corollary 7. For all t ∈ [0,T], we have

|(u− Ū)(t)| + |(v − V̄ )(t)| ≤ CN−1(ln N)2.

Remark 8. According to [7, pg. 186], it is of
interest to generate an accurate approximation
to the rate of reaction u′(t). Given the previous
results, we have that

|(u′−DU)(tj)|≤C|v(tj)−v(tj−1)|+CN−1(ln N)2

≤Ckj‖v′‖(tj−1,tj) +CN
−1(ln N)2.

In the fine mesh

kj‖v′‖(tj−1,tj) ≤ CN
−1 ln N

and in the coarse mesh, where tj > σ,

|v(tj) −v(tj−1)| ≤ C
∣∣ u(tj)
u(tj) + K

−
u(tj−1)

u(tj−1) + K
∣∣

+CN−1 ≤ CN−1,

as ‖u′‖ ≤ C. Hence, we have established the
nodal error bound

|(u′ −DU)(tj)| ≤ CN−1(ln N)2.

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 9 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

0 2 4 6 8 10 12 14 16 18 20

time s

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

U
 a

n
d

 V

epsilon = 2.1362e-05, N= 64

U
V

Figure 1. The numerical solution U, V on the interval [0, 20].

Moreover, as (2.2a) implies that

u′(t) −u′(tj) = (1 −v(tj))(u(tj) −u(t)))
+(u(t) + K−α)(v(t) −v(tj)),

this nodal error bound can be extended to the
global error bound. For all 1 ≤ j ≤ N

|u′(t) −DU(tj)| ≤ CN−1(ln N)2, t ∈ [tj−1, tj].

V. NUMERICAL RESULTS

In this section we present some numerical re-
sults for the numerical method (3.8) applied to the
system (2.2) where the parameters in (2.1) have
been taken to be

k1 = 16847, k−1 = 7, k2 = 12,

s(0) = 2.5 × 10−3, e(0) = 5.4 × 10−8.

This yields the following parameter values for the
non-dimensional system (2.2)

K= 0.4511, α= 0.2849 ε=
e(0)

s(0)
= 1.4 × 2−16.

These values are used in [1] to fit the Henri-
Michaelis-Menten system to experimental data
derived from acetylcholine hydrolysis by acetyl-
cholinesterase. We examine the performance of

the numerical method over an extensive range of
the parameter ε, which is equivalent to varying
the initial concentration e(0) of the enzyme. In
Figure 1 the computed approximations U,V (gen-
erated from the finite difference scheme (3.8)) are
displayed. The plots confirm that v has an initial
layer. Larger values of the parameter ε will result
in less steep gradients appearing in the plot of v.

The global orders of convergence of the finite
difference scheme (3.8) are estimated using the
double-mesh principle [2, Chapter 8, pg. 170].
Note that this principle provides estimates of the
orders of convergence despite the fact that the
exact solution is unknown.

We denote by UN and U2N the computed
solutions on the Shishkin meshes ΩN and Ω2N

respectively. These solutions are used to compute
the maximum two-mesh global differences

D̄Nε := ‖Ū
N − Ū2N‖ΩN∪Ω2N ,

‖g‖ΩN := max
tj∈ΩN

|g(tj)|,

where ŪN (Ū2N ) denotes the linear interpolant of
the discrete solutions UN (U2N ) over the mesh
ΩN (Ω2N ), respectively. The uniform global two-
mesh differences D̄N and their corresponding uni-

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 10 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

Table I
TWO MESH GLOBAL DIFFERENCES D̄Nε , PARAMETER-UNIFORM TWO MESH DIFFERENCES D̄

N AND ORDERS OF
PARAMETER-UNIFORM GLOBAL CONVERGENCE p̄N FOR THE COMPONENT U

ε N=8 N=16 N=32 N=64 N=128 N=256 N=512

20 0.0388 0.0192 0.0095 0.0069 0.0047 0.0029 0.0018

2−2 0.0784 0.0493 0.0154 0.0073 0.0045 0.0026 0.0015

2−4 0.1155 0.0817 0.0300 0.0068 0.0046 0.0036 0.0021

2−6 0.1278 0.0928 0.0361 0.0100 0.0054 0.0042 0.0025

2−8 0.1312 0.0957 0.0377 0.0109 0.0056 0.0044 0.0026

2−10 0.1320 0.0965 0.0381 0.0111 0.0057 0.0045 0.0027

2−12 0.1322 0.0967 0.0382 0.0112 0.0057 0.0045 0.0027

. . . . . . .

. . . . . . .

2−30 0.1323 0.0967 0.0383 0.0112 0.0057 0.0045 0.0027

D̄N 0.1323 0.0967 0.0383 0.0112 0.0057 0.0045 0.0027

p̄N 0.4516 1.3375 1.7734 0.9684 0.3498 0.7491 0.8875

Table II
TWO MESH GLOBAL DIFFERENCES D̄Nε , PARAMETER-UNIFORM TWO MESH DIFFERENCES D̄

N AND ORDERS OF
PARAMETER-UNIFORM GLOBAL CONVERGENCE p̄N FOR THE COMPONENT V

ε N=8 N=16 N=32 N=64 N=128 N=256 N=512

20 0.0428 0.0349 0.0249 0.0163 0.0101 0.0060 0.0035

2−2 0.0573 0.0322 0.0275 0.0207 0.0124 0.0067 0.0034

2−4 0.0484 0.0332 0.0324 0.0260 0.0152 0.0081 0.0042

2−6 0.0454 0.0365 0.0337 0.0272 0.0159 0.0085 0.0045

2−8 0.0461 0.0374 0.0340 0.0275 0.0160 0.0087 0.0045

2−10 0.0463 0.0376 0.0341 0.0276 0.0161 0.0087 0.0045

2−12 0.0463 0.0377 0.0341 0.0276 0.0161 0.0087 0.0046

. . . . . . .

. . . . . . .

2−30 0.0463 0.0377 0.0341 0.0276 0.0161 0.0087 0.0046

D̄N 0.0573 0.0377 0.0341 0.0276 0.0161 0.0087 0.0046

p̄N 0.6039 0.1440 0.3050 0.7781 0.8873 0.9358 0.9663

form orders of convergence p̄N are calculated from

D̄N := max
ε∈S

D̄Nε , p̄
N := log2

(
D̄N

D̄2N

)
,

where S = {20, 2−1, . . . , 2−30}. The maximum
two-mesh global differences D̄Nε for each com-
ponent u,v are, respectively, displayed in Tables I
and II. The uniform two-mesh global differences
D̄N , and their global orders of convergence p̄N

are given in the last two rows of each table. These
numerical results are in line with the asymptotic
error bound established in Theorems 5 and 6. In
the final Table III, the uniform two-mesh global
differences and their global orders of convergence
for the discrete approximations to the rate of
reaction u′(t) are displayed, which again show

parameter-uniform convergence. In all three Tables
we observe that, for N sufficiently large, the global
orders of convergence p̄N are tending towards the
rate associated with the bound N−1(ln N)2.

VI. CONCLUSION

A numerical method is constructed for a singu-
larly perturbed system of two coupled nonlinear
initial value equations. Theoretical error bounds
are established at all time points, which guarantee
that the numerical approximations converge to the
continuous solution, irrespective of the size of the
singular perturbation parameter. Numerical results
support these theoretical error bounds.

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 11 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227


J. J. H. Miller, E. O’Riordan, Robust numerical method for a singularly perturbed problem arising in ...

Table III
TWO MESH GLOBAL DIFFERENCES PARAMETER-UNIFORM TWO MESH DIFFERENCES AND GLOBAL ORDERS OF

PARAMETER-UNIFORM CONVERGENCE FOR THE DISCRETE DERIVATIVE D−U

ε N=8 N=16 N=32 N=64 N=128 N=256 N=512

20 0.0864 0.0841 0.0695 0.0509 0.0340 0.0212 0.0126

2−2 0.1324 0.0706 0.0484 0.0336 0.0248 0.0138 0.0100

2−4 0.1415 0.0736 0.0546 0.0368 0.0268 0.0149 0.0104

2−6 0.1424 0.0747 0.0566 0.0376 0.0274 0.0152 0.0105

2−8 0.1425 0.0754 0.0572 0.0378 0.0276 0.0153 0.0105

2−10 0.1426 0.0756 0.0573 0.0379 0.0276 0.0153 0.0105

2−12 0.1426 0.0757 0.0573 0.0379 0.0276 0.0153 0.0105

. . . . . . .

. . . . . . .

2−30 0.1426 0.0757 0.0573 0.0379 0.0276 0.0153 0.0105

D̄N 0.1426 0.0841 0.0695 0.0509 0.0340 0.0212 0.0126

p̄N 0.7615 0.2742 0.4503 0.5839 0.6812 0.7494 0.7964

REFERENCES

[1] S. Dimitrov, S. Markov, Metabolic rate constants: some
computational aspects, Math. Comput. Simulation, 133,
91–110, (2017).

[2] P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan,
G. I. Shishkin, Robust Computational Methods for
Boundary Layers, Chapman & Hall, New York, 2000.

[3] V. Henri, Lois generales de l’action des diastases, Her-
mann, Paris, 1903.

[4] R. Ishwariya, J. Princy Merlin, J. J. H. Miller, S. Valar-
mathi, A parameter uniform almost first order convergent
numerical method for a non-linear system of singu-
larly perturbed differential equations, BIOMATH 5, 1-8,
(2016).

[5] L. Michaelis, M. L. Menten, Die Kinetik der Invertin-
wirkung, Biochem. Z. 49 333-369, (1913).

[6] J.J.H. Miller, E. O’Riordan and G.I. Shishkin, Fitted
Numerical Methods for Singular Perturbation Problems,
World-Scientific, Singapore (Revised edition), 2012.

[7] J. D. Murray, Mathematical Biology. I An Introduction,
Springer, Third edition, 2001.

[8] L. A. Segel and M. Slemrod, The quasi-steady-state
assumption: a case study in perturbation, SIAM Rev., 31
(3), 446-477, (1989).

[9] A. Zagaris, H. G. Kaper and T. J. Kaper, Fast and slow
dynamics for the computational singular perturbation
method, Multiscale Model. Simul., 2 (4), 613–638,
(2004).

Biomath 9 (2020), 2008227, http://dx.doi.org/10.11145/j.biomath.2020.08.227 Page 12 of 12

http://dx.doi.org/10.11145/j.biomath.2020.08.227

	Introduction
	Continuous problem
	Discrete problem
	Error analysis
	Numerical results
	Conclusion
	References