Uniformly convergent numerical scheme for singularly per- turbed parabolic delay differential equations

This paper deals with numerical treatment of singularly perturbed parabolic differential difference equations having small shifts on the spatial variable. The considered problem contain small perturbation parameter (ε) multiplied on the diffusion term of the equation. For small values of ε the solution of the problem exhibits a boundary layer on the left or right side of the spatial domain depending on the sign of the convective term. The terms involving the shifts are approximated using Taylor’s series approximation. The resulting singularly perturbed parabolic partial differential equation is solved using implicit Euler method in the temporal discretization with exponentially fitted operator finite difference method in the spatial discretization. The uniform stability of the scheme investigated using comparison principle and discrete solution bound by constructing barrier function. Uniform convergence analysis has been carried out. The scheme gives second order convergence for the case ε > N−1 and first order convergence for the case ε ≪ N−1, where N is number of mesh interval. Test examples and numerical results are considered for validating the theoretical analysis of the scheme. 2010 Mathematics Subject Classification: 65M06, 65M12, 65M12.


Introduction
Differential equations with delay terms play an important role in modeling many process in Computational Bio-science, Control Theory, Economics and Engineering [5]. Some of the well known applications are the mathematical modeling of population dynamics and epidemiology [8], physiological kinetics [1], blood cell production [10] and so on. In computational neuroscience, the feasibility of recording single neuron movement induces the development of accurate mathematical models of neuronal variability. In modeling of spiking movement of neuron to any level of exactness, one has to consider special features of each kind of neuron and its input processes [3]. In 1965, Stein [16] developed a mathematical model for the stochastic movement of neuron. After two years, the author [17] generalized his model to handle a distribution of past synaptic potential amplitudes. The Stein's model is once again generalized by Musila and Lansky [11], in the form of singularly perturbed parabolic differential difference equations (SPPDDEs) as to consider the time evolution trajectories of the membrane potential, where µ D and σ are diffusion moments which characterize the influence of dendrites synapses on the cell excitation. a s and λ s are the amplitude and intensity of the excitatory inputs and i s and ω s represent the amplitude and intensity of the inhibitory inputs to the membrane potential. First derivative term represent the exponential decay between two consecutive jumps generated by the input processes. The membrane potential decays to the resting level with time constant τ.
Driving the exact solution of this model problem is very difficult [11]. To simulate the solution, following the numerical techniques becomes mandatory. Since the model problem is in the form of singularly perturbed parabolic differential difference equations, classical numerical methods fail to give a relevant numerical solution. Up to date only a few numerical schemes have been developed for solving this model problem.
In recent years, few scholar's have developed numerical scheme that converges independently of the perturbation parameter for solving SPPDDEs. Ramesh and Kadalbajoo are the first authors, who treat this model problem in time dependent form. In [13], they developed a scheme using backward Euler in time direction and upwind and midpoint upwind FDM on Shishkin mesh in space. Kumar and Kadalbajoo in [9] developed a numerical scheme using backward Euler in time direction with B-spline collocation on Shishkin mesh in space. In series of papers [2][3][4], Bansal and Sharma developed different numerical schemes using non standard FDM for this problem when the delays are large. Gupta et al. in [6] developed scheme using implicit Euler on uniform mesh in time and hybrid finite difference scheme on Shishkin mesh in space and, apply Richardson extrapolation technique for increasing the order of convergence. In [19,20], Woldaregay and Duressa developed uniformly convergent scheme using non standard FDM in space and Runge-Kutta in time. However, the aforementioned methods still lack accuracy which is one of the crucial criteria for the numerical method to be good.
The purpose of this study is to construct more accurate uniformly convergent numerical scheme using exponentially fitted operator finite difference method. Moreover, to discuss the stability and uniform convergence of the scheme. Notation: In this paper N, M denoted the number of mesh intervals in space and time direction respectively. C is the notation for a positive constant independent of c ε , N and M. The norm ∥.∥ is used to denote the maximum norm defined as ∥g∥ = max x,t |g(x, t)|.

Statement of the Problem
The singularly perturbed parabolic differential difference equation of convection diffusion type having reaction term with a retard and advance argument is given as on the domain D = Ω × Λ = (0, 1) × (0, T ] with smooth boundary ∂D =D − D, where T is a fixed positive number, with the initial and interval conditions where ε, 0 < ε ≪ 1 is the singular perturbation parameter and δ, η are the delay and advance parameters satisfying δ, η = o(ε). The functions a(x), α(x), β(x), ω(x), f (x, t), u 0 (x), ϕ(x, t) and ψ(x, t) are assumed to be smooth, bounded and not a function of ε. We assume the coefficients of reaction term α(x), β(x) and ω(x) satisfy for some constant ζ. This condition ensures that the solution of the problem in (2)-(3) exhibits boundary layer in the neighbourhood of D L or D R depending on whether a(x)−δβ(x)+ηω(x) < 0 or > 0 for x ∈Ω. When the shift parameters are zero (i.e. δ = η = 0) the above equation reduces to a singularly perturbed parabolic differential equation (which is one parameter problem), which with small ε exhibits layers depending upon the value of the coefficient function a(x). When a(x) < 0 a regular boundary layer appears in the neighborhood of D L and a(x) > 0 corresponds to existence of a boundary layer near D R . If a(x) changes sign then interior layer (shock layer) will appear on the solution of the problem [6]. The layer is maintained for δ, η 0 but sufficiently small. Note that the above problem in (2)-(3) reduces to the classical case of singularly perturbed parabolic problem when δ, η = 0. The existence of the parameter ε on the diffusion term leads to bad approximation in the computed solution while using standard numerical methods [13]. To avoid these bad approximations, an unacceptably large number of mesh points are required when ε is very small. So, to overcome the drawback associated with standard numerical methods, we developed a numerical scheme using backward Euler in time and exponentially fitted operator FDM in space, which treat the problem very well. The developed scheme is based on the procedure of Rothe's (i.e. first discretizing in time followed by discretization in space).
In the considered case the boundary layer occurs near x = 1. Using compatibility conditions, we deduce that there exist a constant C independent of c ε such that ∀(x, t) ∈D, we have the following conditions that guarantee the existence of a constant C independent of c ε such that ∀(x, t) ∈D for the detail one can refer [14].
Remark 2.1 Note that there does not exist a constant C independent of c ε such that |u(x, t) − ψ(1, t)| ≤ Cx because a boundary layer occurs near x = 1.
The problem obtained by setting c ε = 0 in (5)-(6) is called reduced problem and given as with given initial data It is a first order hyperbolic PDEs, for small values of c ε the solution u(x, t) of (5)-(6) is very close to the solution u 0 (x, t) of (7)-(8).

Properties of analytical solution
Since u 0 (x) is sufficiently smooth and using the property of norm, we prove the following lemma. (5)- (6). Then we obtain the bound

Lemma 2.3 Let u(x, t) be the solution of the continuous problem in
where ζ > 0 is lower bound of q(x) and |u(x, t)| ∂D is the restriction of |u(x, t)| on ∂D.

Proof. By defining barrier functions
and applying the maximum principle we obtain the required bound.

Lemma 2.4 The bound on the derivative of the solution u(x, t) of the problem in (5)-(6) with respect to x and t is given by
where p * is lower bound of p(x).
Proof. Assume x * ∈ [0, 1] be such that z j+1 (x * ) = min x∈Ω z j+1 (x) < 0. From the above assumption, it is clear that x * {0, 1} implies that x * ∈ (0, 1). Applying the property in calculus, we have d dx z j+1 (x * ) = 0 and d 2 Next, we analyze the discretization error, let denote the local truncation error by e j+1 := Lemma 3.2 The local truncation error estimate in the temporal direction is given by Proof. Using the Taylor series expansion to u(x, t j+1 ), we obtain Substituting (15) into (5) gives Since error satisfies the DEs, the local truncation error satisfies Using the semi-discrete maximum principle gives Next, we need to give the bound for the global truncation error of the semi-discretization. Let denote T E j+1 be the global error estimate up to the ( j + 1)th time step.

Lemma 3.3 The global error estimate at t j+1 is given by
Proof. Using the results in local error estimate up to the ( j + 1)th time step in Lemma 3.2, the global error is given as where C is constant independent of c ε and ∆t.
Next we set a bound for the derivatives of solution of the problem in (12).  (12)- (13) satisfies the bound Proof. Direct from Lemma 2.4.

Spatial Discretization
Generally, there are two major techniques for designing numerical methods which gives small truncation errors in the boundary layer. The first approach is the class of layer adapted mesh methods which uses fine mesh in the layer region and coarse mesh in outer layer region. The second approach is the exponentially fitted operator methods, in which it uses uniform mesh and an exponentially fitting factor is determined and induced on the diffusion term of the equation for stabilizing the term containing the singular perturbation parameter. In this approach the difference schemes reflect the qualitative behaviour of the solution inside the boundary layer region.
In this paper, we use exponentially fitted operator finite difference method (FOFDM) which helps us to hinder the influence of the singular perturbation parameter(ε) as ε → 0 in the boundary layer region. For each j = 0, 1, 2, ..., M − 1 we have the boundary value problems of the form with boundary conditions Using the theory applied in asymptotic method we develop exponentially fitted numerical scheme to solve the singularly perturbed BVPs in (20)-(21). In the considered case the boundary layer is in the right side of the domain, so the zeroth order asymptotic solution of (20)-(21) is given as where U j+1 0 is the solution of the reduced problem [12]. The spatial domainΩ = [0, 1] is discretized into N equal number of subintervals each of length h = N −1 . Let 0 = x 0 , x N = 1 and x i = ih, i = 0, 1, 2, ..., N be the mesh points.
Considering h is small enough, the discretized form of (22) becomes Similarly, we write where ρ = h/c ε , h = 1/N. To handle the effect of the perturbation parameter artificial viscosity (exponentially fitting factor σ(ρ)) is multiplied on the term containing the perturbation parameter as Next, on a uniform points {x i } N i=0 with h = x i+1 − x i . For any mesh function Z i , we denote the difference operators as Applying the central finite difference on (25) gives Multiplying (27) by h and considering h → 0 and truncating the term Substituting the results in (23) and (24) into (28) and simplifying, the required exponential fitting factor is obtained as Hence, the proposed finite difference scheme becomes where

Stability and Uniform Convergence Analysis
The proposed discrete scheme satisfies the discrete comparison principle.  [7].
The result of this lemma guarantee the existence and uniqueness of the discrete solution.
Proof. The proof is a simple computation, enables one to give a bound, that is uniform in c ε for the norm of the inverse of L h,∆t .

Lemma 3.7 [Uniform stability estimate.] The solution U j+1
i of the discrete scheme in (30) satisfy the following bound.
Proof. Consider two barrier functions of the form π ± i, j+1 = 1 Hence, using the discrete comparison principle gives This result guarantee the uniform stability of the discrete scheme. Next, we consider the semi-discrete problem in (12)-(13) and the fully discrete scheme in (30) to estimate the truncation error of the spatial direction discretization.

Using the bounds for the derivatives of the solution in Lemma 3.4 gives
which gives the required bound. Proof. Refer on [15]. Using Lemma 3.9 into (36) gives Hence, using the bound in Lemma 3.8 gives Theorem 3.2 Let u(x i , t j+1 ) and U j+1 i be solution of (5)- (6) and (30) respectively on discretized domain. Then for sufficiently large N, the following uniform error estimate holds: Proof. The combination of temporal and spatial error bound gives the required result.

Numerical Examples and Results
In this section, we illustrate the proposed scheme using two numerical examples of the form given in (2)-(3). The exact solution of these two examples are not known. We investigate the theoretical results in this paper by performing experiments using the proposed scheme.

Example 4.1 Consider the problem
with the initial condition u 0 (x) = 0, x ∈ [0, 1] and the interval conditions ϕ( into the mesh points. We calculate the maximum absolute error as

Example 4.2 Consider the problem
and the ε-uniform (parameter uniform) error estimate as We calculate the rate of convergence of the proposed scheme using ε,δ,η and the εuniform rate of convergence using 1.1208e+00 1.1071e+00 1.0170e+00 8.2794e-01 7.3438e-01 6.9902e-01 2 −10 1.1316e+00 1.1463e+00 1.1470e+00 1.1219e+00 1.0283e+00 8.4274e-01 2 −12 1.1323e+00 1.1489e+00 1.1568e+00 1.1589e+00 1.1533e+00 1.1254e+00 2 −14 1.1324e+00 1.1490e+00 1.1574e+00 1.1614e+00 1.1629e+00 1.1619e+00 2 −16 1.1324e+00 1.1490e+00 1.1575e+00 1.1615e+00 1.1635e+00 1.1643e+00 2 −18 1.1324e+00 1.1490e+00 1.1575e+00 1.1615e+00 1.1635e+00 1.1644e+00 2 −20 1.1324e+00 1.1490e+00 1.1575e+00 1.1615e+00 1.  Figure 1 and 2, as ε goes very small. In Figure 3, effect of the delay parameter on the behaviour of the solution is shown. When the magnitude of delay parameter grows, the thickness of the boundary layer decreases as seen on Figure 3 (a) and (b). The maximum absolute error, rate of convergence, ε-uniform error and ε-uniform rate of convergence of Example 4.1 is given in Table 1-3 and Example 4.2's in Table 5-7. The result in Table 1 shows, as the perturbation parameter ε 2 > N −1 the scheme gives good approximation to the exact solution, whereas as ε goes small the numerical scheme without the exponential fitting factor diverges. But, as we observe the results in Table 2, the exponentially fitted scheme converges independent of the perturbation parameter as ε → 0. Similarly, we observe the same scenario in the    Table 5 and 6. In Table 3 and 7 we observe, the maximum absolute error and rate of convergence of the scheme for different values of delay and advance parameter by keeping ε = 2 −2 > N −1 . As one observes in these tables the scheme    4, the Log-Log plot of the maximum absolute error verses N are given for singular perturbation parameter ranging from ε = 2 −4 to 2 −20 . In this figure the graphs are parallel and overlapped as ε goes small, this indicate that the proposed scheme converges independent of the values of perturbation parameter. This is one of the main result claimed to be shown in this paper. As one observe in Table 2 and 6, for each N and M as the perturbation parameter (ε) → 0, the maximum absolute error is stable and uniform. We compared the ε-uniform error and ε-uniform rate of convergence of the proposed scheme with some published papers in Table 4 and 8, as one observes the proposed scheme gives more accurate result than schemes in [9], [13] and [20].

Conclusion
Singularly perturbed parabolic differential difference equations exhibiting boundary layer is considered. The considered problem contains small perturbation parameter multiplied to the highest order derivative term of the equation; and small delay and advance parameters on the non derivative terms of the spatial variable. Exponentially fitted operator numerical scheme is proposed for solving the problem. First, the equation is approximated by equivalent singularly perturbed parabolic PDEs using Taylor's series approximation for the terms with the delay and advance parameters. Inducing exponential fitting factor for term with perturbation parameter and determine its value. Numerical scheme is developed using implicit Euler in temporal discretization with central FDM for spatial discretization. The uniform stability and uniform convergence of the scheme is established. It is shown that the scheme is accurate and converges uniformly with order of convergence O(N −1 + ∆t).