A finite difference method on layer-adapted mesh for singularly perturbed delay dif- ferential equations

The purpose of this paper is to present a uniform finite difference method for numerical solution of a initial value problem for semilinear second order singularly perturbed delay differential equation. A numerical method is constructed for this problem which involves appropriate piecewise-uniform Shishkin mesh on each time subinterval. The method is shown to uniformly convergent with respect to the perturbation parameter. A numerical experiment illustrate in practice the result of convergence proved theoretically.

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 where δ 0 , c 0 and c 1 are nonnegative constants. Then it holds that Proof.The semilinear equation (1) can be written in the form where We rewrite (11) in the form where From (12), we have the following relation for u (t) Integrating (13) from 0 to t, we have By the change of integral bounds, we get Applying Lemma 2.1, we obtain (7). Due to F ∞ ≤ C, it now follows from (13) that where From (14) we have the following relation for u (t) and it is easy to see that It follows from (1.1) that Using the estimates (8), (16) and (17) in (21) for k = 0, we have i.e., the inequality (9) is valid for p = 1.
Next, from the last inequality for t = r, we have u (r) ≤ C and thereby from (15) it follows that i.e. the inequality (9) is also valid for p = 2. Since 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ω N 0 be any non-uniform mesh onĪ: which contains by N mesh point at each subinterval I p (1 ≤ p ≤ m) : and consequently To simplify the notation we set w i = w(t i ) for any function w(t), and moreover y i denotes an approximation of u(t) at t i . For any mesh function {w i } defined onω N 0 , we use For the difference approximation of (1), we are using the following identity where exponential basis functions We note that functions ψ 1i (t) and ψ 2i (t) are the solutions of the following problems, respectively: Using interpolating quadrature rules with the weight and remainder terms in integral form on subintervals [t i−1 ,t i+1 ], consistent with [23,24], after a simple calculation, we have the following relation: where the factor and remainder term To define an approximation for the boundary condition (3), we proceed our discretization process by where the exponential basis function We note that functions ϕ 0 (t) is the solution of the following problem: Whence, as similar above we can write the following difference relation: where the coefficient and the remainder term Neglecting R i and r (0) in (19) and (22), we propose the following difference scheme for approximation (1)- (2): 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 [r p−1 , σ p ] and [σ p , r p ] into N/2 equidistant subintervals, where the transition point σ p , which separates the fine and coarse portions of the mesh is obtained by where β 1 ≥ 1 and β p > 1 (2 ≤ p ≤ m) are some constants. Hence, if h (1) p and h (2) p denote the stepsizes in [r p−1 , σ p ] and [σ p , r p ] respectively, we have In the rest of the paper we only consider this type mesh.

Convergence Analysis
We now estimate the approximate error z i = y i − u i , 0 ≤ i ≤ N 0 , which satisfies the discrete problem where the truncation error R i and r (0) are given by (21) and (24).
where θ * = ρα 2 coth(αρ/2), γ = 4α −1 exp(4α −1 (b * + c * )) Proof. Let z t,i = v i . Then the relation (27) can be rewritten as Solving the first order difference equation with respect to v i , we get where Then, since and taking into consideration (5), from (30) after some manipulations we have the following inequality Therefore 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 ∈ C 1 (Ī), ψ ∈ C 1 (Ī 0 ). Then for the truncation errors R i and r (0) , the following estimates hold We now consider the case σ p = r p−1 + α −1 β p ε ln N and estimate R i on [r p−1 , σ p ] and [σ p , r p ] separately. In the layer region [r p−1 , σ p ], the inequality (32) reduces to Hence Next we estimate R i for (p − 1/2)N + 1 ≤ i ≤ pN. In this case, recalling that p , we obtain from (32) and this implies that 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 Combining the two previous lemmas gives us the following main convergence result. Theorem 4.1. Let a ∈ C 1 (Ī), ϕ ∈ C 1 (I 0 ) . The continuously differentiable function f (t, u, v) satisfies the conditions (1.4) and the derivative ∂ ∂t f (t, u, v) is bounded for 0 ≤ t ≤ T and |u|, |v| ≤ C. Then the following estimate holds

Numerical Results
In this section, a simple numerical example is devised to verify the validity for the proposed method. Consider the test problem with 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 N,p 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.  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.