UNIFORMLY CONVERGENT NUMERICAL METHOD FOR SINGULARLY PERTURBED DELAY PARABOLIC DIFFERENTIAL EQUATIONS ARISING IN COMPUTATIONAL NEUROSCIENCE

The motive of this work is to develop ε-uniform numerical method for solving singularly perturbed parabolic delay differential equation with small delay. To approximate the term with the delay, Taylor series expansion is used. The resulting singularly perturbed parabolic differential equation is solved by using non-standard finite difference method in spatial direction and implicit Runge-Kutta method for the resulting system of IVPs in temporal direction. Theoretically the developed method is shown to be accurate of order O(N−1 + (∆t)2) by preserving ε-uniform convergence. Two numerical examples are considered to investigate εuniform convergence of the proposed scheme and the result obtained agreed with the theoretical one.


Introduction
The class of differential difference equations which have characteristics of delay and singularly perturbed behavior is known as singularly perturbed differential difference equations or singularly perturbed delay differential equations (SPDDEs). SPDDEs plays an important role in modeling real life phenomena in Bioscience, Control Theory, Economics and Engineering [6]. A few application areas are the mathematical modeling of: epidemiology and population dynamics [11], physiological kinetics [1], production of blood cell [16]. Singularly perturbed differential equations relate an unknown function to its derivatives evaluated at the same instance. Whereas, singularly perturbed delay differential equations model physical problems for which the evaluation depend on the present state of the system and also on the past history. In general, when the perturbation parameter approaches to zero the smoothness of the solution of the SPDDEs deteriorates and it forms a boundary layer(s) [19]. The time dependent problem i.e. singularly perturbed parabolic delay differential equations (SPPDDEs) have applications in the study of modeling of neuronal variability [23].
Various mathematical models are given in books [25,26], for the determination of the behavior of a neuron to random synaptic inputs. In 1991, Musila and Lansky [18] gives generalization for the Stein's model [23] and proposed a model to treat the time evolution trajectories of the membrane potential in terms of SPPDDEs as where µ D and σ are diffusion moments of Wiener process characterizing the influence of dendritic synapses on the cell excitability. The excitatory input contributes to the membrane potential by an amplitude a s with intensity λ s and similarly the inhibitory input contributes by an amplitude i s with intensity ω s . The first derivative term is because of the exponential decay between two consecutive jumps caused by the input processes. The membrane potential decays exponentially to the resting level with a membrane time constant τ . This model makes available time evolution of the trajectories of the membrane potential. The model in (1.1) is a singularly perturbed parabolic differential difference equation, one can hardly derive its exact solution.
Hence, to simulate this model, one has to land to a very suitable numerical methods.
In the last few years, scholar's have devoted for the development of numerical solution of this problem. In addition to that, authors try to show the influence of the delay parameter on the behavior of the solution. In research articles [2-4, 12, 13, 19, 20] authors develop different numerical schemes to study a classes of SPPDDEs and discussed the effect of shifts on the solution.
In this paper, we construct and analyze a non-standard finite difference scheme which utilizes on uniform mesh. The proposed scheme uses the procedures of method of line, which consist of non-standard finite difference operator for the spatial discretization and Runge-Kutta method for the time discretization. For the theoretical analysis, the global error is decomposed into two parts: the first is due to the spatial discretization and the second is due to the temporal discretization of the semi-discrete problem obtained after the spatial discretization.
The main contribution of this work is, to develop ε-uniform numerical scheme for the SPPDDEs containing small delay on the spatial variable. In the proposed method, it is not required to have any restriction on the mesh generation.
This paper is organized as follows. In Section 1 a brief introduction about the problem is given, in Section 2 definition of the problem and the behavior of its analytical solution is given. In Section 3, discretizing the spatial domain and techniques of non-standard finite difference is discussed, and the ε-uniform convergence of the semidiscrete problem is proved. Next, Runge-Kutta method used for the system of IVPs resulted from spatial discretization and discuss the convergence of the discrete scheme.
In Section 4, numerical examples and results are given to validate the theoretical analysis and finally, in Section 5, the conclusion of the work done is presented.
Notations. Through out this paper N, M denoted for the number of mesh elements in space and time direction respectively. C (in some case indexed) is denoted for positive constant independent of perturbation parameter and N . The norm ∥·∥ Ω N x ×Ω M t is used to denote discrete maximum norm.

Problem Formulation
A class of singularly perturbed parabolic differential difference equation with delay on the spatial variable on domain D = Ω x × Ω t is given by where (x, t) ∈ D = Ω x × Ω t = (0, 1) × (0, T ] for fixed positive number T, ε is a singular perturbation parameter with 0 < ε ≪ 1 and δ is delay parameter assumed to be sufficiently small as order of o(ε). We assume, the functions a(x), α(x), β(x) and f (x, t), u 0 (x), ϕ(x, t) and ψ(x, t) are sufficiently smooth, bounded and independent of ε. The coefficients of reaction term β(x) and delay term α(x) are assumed to satisfy α(x) + β(x) ≥ θ > 0, for all x ∈Ω x , for some positive constant θ. This condition ensure that the solution of (2.1) exhibits boundary layer in the neighborhood of When the shift parameter is zero (i.e., δ = 0) the above equation reduces to a singularly perturbed parabolic differential equation, which exhibits layer(s) depending on the value of a(x). If a(x) < 0 a boundary layer appears in the neighborhood of x = 0, in case a(x) > 0 corresponds to existence of a boundary layer near x = 1. For the case a(x) change sign interior layer will appear on the solution of the problem [9]. The layer is maintained for δ ̸ = 0 but sufficiently small.
When the shift parameter δ < ε, the use of Taylor's series expansion for the term containing shift is valid [24].
2.1. Estimate for the delay term. From the assumption δ < ε, by using Taylor series approximation for u(x − δ, t), we obtain Now, substituting this approximations into (2.1), we obtain For small δ, equations in (2.1) and (2.3) are asymptotically equivalent, because the difference between the two equations is order of O(δ 3 ). Now we assume again 0 We also assumed that p(x) = a(x)−δα(x) ≥ p * > 0, which guarantee the occurrence of boundary layer in the neighborhood of x = 1. The other case p(x) = a(x) − δα(x) ≤ p * < 0, imply the occurrence of the boundary layer in the neighborhood of x = 0 and can be treated in similar manner.
We impose the compatibility conditions so that the data matches at the two corners (0, 0) and (1, 0). The considered case is when boundary layer occurs near x = 1. Now, using compatibility conditions in (2.4) and (2.5), we have the conditions that guarantee the existence of a constant C independent of c ε such that for all (x, t) ∈D for the detail of this one can see [21, page 105]. By setting c ε = 0 in (2.3) we obtain the reduced problem as: This is a first order hyperbolic PDE with initial data given along the two sides t = 0 and x = 0 of the domainD. For small values of c ε the solution u(x, t) of the problem in (2.3) will be very close to u 0 (x, t). To obtain error bounds on the solution of the difference scheme, we assume that the solution of the reduced problem in (2.6) is sufficiently smooth.

Properties of continuous solution.
To show a boundedness of the solutions u(x, t) of (2.3), we assume the initial condition to be zero. Since u 0 (x) is assumed sufficiently smooth and using the property of norm, we prove the following lemma.

and implies that
Lv(x * , t * ) < 0, which contradict to the assumption made above. Here, we have

Lemma 2.3 (Stability estimate). Let u(x, t) be the solution of the continuous problem in (2.3). Then we have the bound
Proof. We define two barrier functions ϑ ± as At the initial value At the boundary points Hence, the proof is completed.

Formulation of Numerical Scheme
3.1. Discretization in spatial direction. The theoretical basis of non-standard discrete modeling method is based on the development of exact finite difference method. In [17], Micken's presented techniques and rules for developing non-standard FDMs for different problem types. In Mickens's rules, to develop a discrete scheme, denominator function for the discrete derivatives must be expressed in terms of more complicated functions of step sizes than those used in the standard procedures. These complicated functions constitutes a general property of the schemes, which is useful while designing reliable schemes for such problems.
For the problem of the form in (2.3), in order to construct exact finite difference scheme we follow the procedures used in [3]. Consider the constant coefficient subequations given in (3.1) and (3.2) by ignoring the time variable as where p(x) ≥ p * and q(x) ≥ θ. Thus the problem (3.1) has two independent solutions namely exp(λ 1 x) and exp(λ 2 x) with We discretize the spatial domain [0, 1], using uniform mesh length ∆x = h such that where N is the number of mesh points in spatial direction. We denote the approximate solution of u(x) at x i 's by U i . Now our objective is to calculate a difference equation which has the same general solution as the problem (3.1) has at the grid point x i given by Simplifying (3.4), we obtain that is an exact difference scheme for (3.1). After doing the arithmetic manipulation and rearrangement on (3.5) we obtain The denominator function becomes γ 2 = hcε p * exp hp * cε − 1 . Adopting this denominator function for the variable coefficient problem, we write it as where γ 2 i is a function of c ε , p i and h, where p i is denoted for p(x i ). By using (3.7) in to the main semi-descrete scheme, we obtain Let U i (t) is denoted for the approximation of u(x i , t). By using the non-standard finite difference approximation. At this stage the problem in (2.3) reduces to semi-discrete form as The system of IVPs in (3.8) can be written in compact form as where A is a tridiagonal matrix of size N − 1 × N − 1 and U i (t) and F i (t) are column vectors with N − 1 entries. The entries of A and F are given respectively as Now we need to show the semi-discrete operator L h satisfies the maximum principle and the uniform stability estimate.

Lemma 3.2. The solution U i (t) of the semi-discrete problem in
. At the boundary points From Lemma 3.1, we obtain ϑ ± i ≥ 0, for all (x i , t) ∈Ω N x × Ω t . This complete the proof. □ 3.2. Convergence Analysis. Now, let us analyze the convergence of the spatial discretization. We proved above the semi-discrete operator L h satisfy the maximum principle and the uniform stability estimate. Note that U i (t) is denoted for the spatial semi-discretization approximate solution to the exact solution u(x, t) at x = x i , i = 0, 1, . . . , N . Let define the forward and backward finite differences in space as: respectively, and the second order finite difference operator as Proof. We consider the truncation error in spatial discretization as: The bound c ε h 2 γ 2 i − 1 ≤ Ch used in above expression is based on the behavior of the denominator function (γ 2 i ) in non-standard finite difference. To illustrate the bound given there, let us define σ =: By simplifying and writing explicitly we obtain and we obtain the limit is bounded as lim σ→0 R(σ) = 1 2 , lim σ→∞ R(σ) = 0. Hence, for all σ ∈ (0, ∞), we have R(σ) ≤ C. So, the error estimate in the spatial discretization is bounded as Proof. Let us consider the partition of the interval [0, 1], with points x j = jh, h = 1/N, j = 0, 1, . . . , N . By using

From (3.12) and boundedness of derivatives of solution in Lemma 2.4, we obtain
Then, by the repeated application of L'Hospital's rule, we have Hence, the proof is completed. □ (3.14) sup Proof. Results from boundedness of solution, Lemma 3.1 and Theorem 3.1 gives the required estimates.
with the initial condition U ( . Now, we solve the system of IVPs in (3.15) using implicit Runge-Kutta method of order 2 and 3 given in [7].

Lemma 3.4. From the approximation method in temporal direction, the global error estimates in this direction is given by
where ∥E j+1 ∥ is the global error in the temporal direction at (j + 1) th time level.
Proof. Using the local error estimate upto jth time step, we obtain the global error estimate at (j + 1)th time step. Hence, Here, since C and ∆t are independent of the perturbation parameter ε. By taking the suprimum for all ε ∈ (0, 1] we obtain This shows that the discretization in time direction is consistent and global error is bounded, with the error bound O(∆t) 2 . Now, we use (3.14) and (3.18) to prove the ε-uniform convergence of the fully discrete scheme as sup Hence, we obtain the required bound as (3.19) shows the ε-uniform convergence of the scheme with order of convergence: first order in spatial direction and second order in temporal direction.

Numerical Results and Discussion
To validate the established theoretical results, we perform some experiments using the proposed numerical scheme on the problem given in (2.1). We consider two numerical examples to verify the ε-convergence of the proposed scheme. Example 4.2. From [9] by setting η = 0 with T = 2 subject to the initial condition u(x, 0) = 0, x ∈ [0, 1], and interval-boundary Exact solution is not available for these two problems, therefore maximum nodal errors are calculated by using the double mesh technique given in [9] as The rate of convergence of the method is given by ε,δ log 2 and the ε-uniform convergence is calculated by decreases as N increases for each value of ε. in addition to that, the maximum point-wise error is stable as ε → 0 for each N , M . Using the results in these two examples we confirm the proposed numerical method is more accurate, stable and ε-uniform convergent with rate of convergence one. Numerical results in tables and figure confirm the parameter-uniformness of the proposed scheme. The results in this method are better than that obtained in [20] and [19].

Conclusions
A numerical method is developed for solving a singularly perturbed parabolic delay differential equation whose solution exhibit a boundary layer behavior. The developed method is based on method of line that constitute the non-standard finite difference for the spatial discretization and the Runge-Kutta order 2 and 3 implicit method in the temporal direction for the system of IVP resulting from the spatial discretization. Stability and convergence analysis of the proposed scheme is shown. The applicability