Novel Numerical Scheme for Singularly Perturbed Time Delay Convection-Diffusion Equation

This paper deals with numerical treatment of singularly perturbed parabolic di ﬀ erential equations having large time delay. The highest order derivative term in the equation is multiplied by a perturbation parameter ε , taking arbitrary value in the interval ð 0, 1 (cid:2) . For small values of ε , solution of the problem exhibits an exponential boundary layer on the right side of the spatial domain. The properties and bounds of the solution and its derivatives are discussed. The considered singularly perturbed time delay problem is solved using the Crank-Nicolson method in temporal discretization and exponentially ﬁ tted operator ﬁ nite di ﬀ erence method in spatial discretization. The stability of the scheme is investigated and analysed using comparison principle and solution bound. The uniform convergence of the scheme is discussed and proven. The formulated scheme converges uniformly with linear order of convergence. The theoretical analysis of the scheme is validated by considering numerical test examples for di ﬀ erent values of ε .


Introduction
A large number of mathematical models appear in different areas of science and engineering such as in control theory, epidemiology, and laser optics that take into account not only the present state of a physical system but also its past history [1]. These models are described by certain classes of functional differential equations often called delay differential equations. Examples of delays include the time taken for a signal to travel to the controlled object, driver reaction time, the time for the body to produce red blood cells, and cell division time in the dynamics of viral exhaustion or persistence. In the life sciences, delays are often introduced to account for hidden variables and processes which, although not well understood, are known to cause a time lag [1]. Time delays are natural components of the dynamic processes of biology, ecology, physiology, economics, epidemiology, and mechanics [2], and to ignore them is to ignore reality [1].
Singularly perturbed delay differential equations are differential equations in which its highest order derivative term is multiplied by small perturbation parameter ε and having at least one delay term. A simplified mathematical model for the control system of a furnace used to process metal sheets is represented by singularly perturbed parabolic delay differential equation of the form [3] where u is the temperature distribution in a metal sheet, moving at an instantaneous material strip velocity v and heated by distributed temperature source specified by the function f ; both v and f are dynamically adapted by a controlling device monitoring the current temperature distribution. The finite speed of the controller, however, introduces a fixed delay of length τ. When τ = 0, this problem becomes a thermal problem without time delay.
The presence of the singular perturbation parameter ε on the highest order derivative term leads to occurrences of oscillations in the computed solution while using central FDM, Galerkin FEM, and collocation methods on a uniform mesh [4]. While using these methods to avoid the oscillations, an unacceptably large number of mesh points are required when ε is very small. This is not practical due to the limited size of computer memory and rounding off error. Therefore, to overcome this drawback associated with classical numerical methods, authors use the fitted mesh technique which have a fine mesh in the boundary layer region or fitted operator technique. Recently, numerical treatment of singularly perturbed parabolic PDEs has gotten great attention, and different authors in [5][6][7][8] have developed uniformly convergent numerical schemes (scheme that converges independent of the influence of the singular perturbation parameter). Currently, authors in [9][10][11][12][13] have developed numerical schemes for solving singularly perturbed parabolic time delay reaction-diffusion equations. In papers [14,15], the authors considered a singularly perturbed parabolic convection diffusion equation with time delay and degenerate coefficient. They treat the problems using fitted mesh techniques. Different authors in [16][17][18][19][20][21][22][23] developed numerical schemes using fitted mesh techniques for treating singularly perturbed parabolic time delay convection-diffusion equations. In [24], Podila and Kumar used nonstandard FDM for treating singularly perturbed parabolic time delay convection-diffusion equations. Their scheme gives linear order of convergence.
The main motive of this paper is to formulate an accurate and uniformly convergent numerical scheme for the singularly perturbed parabolic time delay convection-diffusionreaction equation using exponentially fitted operator FDM and to establish the stability and uniform convergence of the scheme. The proposed scheme uses the procedure of Roth (i.e., first discretizing in temporal direction followed by discretization in spatial direction) using the Crank-Nicolson method in temporal direction and exponentially fitted operator FDM on spatial direction. In this method, it is not required to have any restriction on the mesh generation. Notation 1. Throughout this paper, N and M are denoted for the number of mesh intervals in space and time direction discretization, respectively; C is denoted for positive constant independent of ε, N, and M. The norm k:k denotes the maximum norm defined as kgk = max x,t jgðx, tÞj.

Continuous Problem
Consider a class of second order singularly perturbed parabolic convection-diffusion-reaction equations having a term with large time delay, on the domain D = Ω x × Ω t = ð0, 1Þ × ð0, T and on smooth boundary ∂D = D l ∪ D b ∪ D r is given by with the boundary conditions and the interval condition where ε ∈ ð0, 1 is a singular perturbation parameter and τ > 0 is the delay parameter. The terminal time T is assumed to satisfy the condition T = kτ for some positive integer k. The functions aðxÞ, bðx, tÞ, cðx, tÞ, and f ðx, tÞ and ϕ l ðtÞ, ϕ b ðx , tÞ, and ϕ r ðtÞ are assumed to be sufficiently smooth and bounded and satisfy This condition ensures that the solution of the problem in (2)-(4) exhibits a boundary layer of thickness OðεÞ on the right side of the spatial domain.
Our objective in this paper is to formulate an accurate and parameter uniformly convergent numerical scheme and to discuss the uniform stability and the parameter uniform convergence of the scheme for the considered problem in (2)-(4).

Bounds on the Solution and Its
Derivatives. The existence and uniqueness of the solution of (2)-(4) can be established by assuming that the data is Holder continuous and imposing appropriate compatibility conditions at the corner points, using the assumptions of sufficiently smoothness ofϕ l ðtÞ, ϕ b ðx, tÞ, andϕ r ðtÞ.
The required compatibility condition at the corner points and the delay term are so that the data matches at the corner points. For aðxÞ, bðx, tÞ, cðx, tÞ, and f ðx, tÞ to be continuous on domain D, then (2)-(4) has the unique solution uðx, tÞ ∈ C 2 ðDÞ ∩ C 0 ð DÞ. In particular, when the compatibility conditions are not satisfied, a unique classical solution still exists but is not differentiable on all of ∂D.
The reduced problem corresponding to the singularly perturbed parabolic delay PDE (2)-(4) is given as Advances in Mathematical Physics with the boundary conditions and the interval condition It is in the form of hyperbolic delay PDEs. The solution uðx, tÞ of (2)-(4) becomes very close to the solution u 0 ðx, tÞ of (8)-(10) as ε ⟶ 0.
Proof. The result follows from the compatibility condition. See the detailed proof in [16].
Lemma 4. Let uðx, tÞ be the solution of the continuous problem in (2)-(4). Then, we obtain the bound where ζ is the lower bound of bðx, tÞ.

Lemma 5.
The bound on the derivative of the solution uðx, tÞ of the problem in (2)-(4) with respect to x and t is given by where α is lower bound of aðxÞ.

Numerical Scheme
In general, for singularly perturbed problems, there are two strategies for designing numerical methods which have a small error in the boundary layer region [25]. The first approach is the class of fitted mesh methods which uses fine mesh in the boundary layer region and coarse mesh in outer layer region. The stability and convergence analyses of this approach are well developed. The second approach is the fitted operator methods in which it uses uniform mesh and an exponentially fitting factor 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.
Next, let us analyse the truncation error for the temporal discretization made above.
The local truncation error in the temporal direction is given by Proof. Using Taylor's series approximation for uðx, t j Þ and uðx, t j+1 Þ centering at t j+1/2 , we obtain From the approximation in (20), we obtain Using the approximation into (2), we obtain where Since the error e j+1 ðxÞ ≕ uðx, t j+1 Þ − U j+1 ðxÞ satisfies the semidiscrete differential operator Hence, by applying the maximum principle, we obtain Next, we need to show the bound for the global error of the temporal discretization. Let us denote TE j+1 as the global error up to the ðj + 1Þth time step. 4 Advances in Mathematical Physics Lemma 8. The global error at t j+1 is given by Proof. Using the local error up to the ðj + 1Þth time step given in the above lemma, we obtain the global error at the ðj + 1Þth time step as where C is a constant independent of ε and Δt. See the detailed proof in [26].
Next, we set a bound for the derivatives of solution of (15)- (16).
Proof. See the proof in [5].
3.2. Spatial Discretization. The spatial domain ½0, 1 is discretized into N equal number of subintervals, each of length h. Let 0 = x 0 , x N = 1, and x i = ih, i = 0, 1, 2, ⋯, N, be the mesh points. For spatial discretization, we apply an exponentially fitted operator finite difference method which helps us to hinder the influence of the singular perturbation parameter. First, let us find the exponential fitting factor for anonymous BVPs and then apply discretization in the spatial direction.
3.2.1. Computing the Exponential Fitting Factor. To develop the numerical method for (15)-(16), we use the technique designed in the theory of asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer occurs on the right side of the domain. From the theory of singular perturbation in [27], the zero order asymptotic solution of the singularly perturbed boundary value problems of the form is given by Using Taylor's series expansion for aðxÞ and bðxÞ about x = 1 and restriction to their first terms and simplifying gives where u 0 is the solution of the reduced problems (obtained by setting ε = 0) in (29). Considering h is reasonably small and evaluating the result in (31) at x i gives Similarly, we have where ρ = h/ε. Consider a uniform grid Ω N x = fx i g N i=0 and denote h = x i+1 − x i . For any mesh function v i , define the following difference operators: To handle the effect of the perturbation parameter, we multiply artificial viscosity (exponentially fitting factor σðρÞ) on the diffusive part of the problem. Introducing an exponentially fitting factor σðρÞ in (29) and applying the central finite difference scheme gives From (32), we have Substituting (37) and (33) into (36) and simplifying, the exponential fitting factor is obtained as 3.2.2. The Discrete Scheme. Using the central finite difference method for the spatial discretization of (15)- (16) and applying the exponential fitting factor in (38), for i = 1, 2, ⋯, N − 1, the fully discrete scheme becomes where In explicit form, the scheme is rewritten as where 3.3. Stability and Uniform Convergence Analysis. First, we need to prove the discrete comparison principle for the discrete scheme in (39).

Lemma 10. Discrete comparison principle. There exist a comparison function
Proof. The matrix associated with operator ð1 + ðΔt/2ÞL Δt,h Þ U i,j+1 is of size ðN + 1Þ × ðN + 1Þ with its entries for i = 1, 2, ⋯, N − 1 are So, the coefficient matrix satisfies the property of M matrix. So, the inverse matrix exists and it is nonnegative. This guarantees the existence and uniqueness of the discrete solution. See the detailed proof in [28].
Proof. The proof is a simple computation, enables one to give a bound, that is uniform in ε for the norm of the inverse of 6 Advances in Mathematical Physics where bðx i , t j+1/2 Þ ≥ ζ > 0.
Proof. On the discrete domain fx i g N 0 , for the interior grid points, we have Since x 1 = h, 1 − x N−1 = h: Repeatedly using L'Hospital's rule gives This completes the proof.
Theorem 16. The numerical solution U i,j+1 of the problem in (39) satisfies the following uniform error bound Proof. Using Lemma 15 into (56) gives

Advances in Mathematical Physics
Hence, using the result in Lemma 13, we obtain Using the supremum over all ε ∈ ð0, 1, we obtain Remark 17. For the case ε > h, the developed scheme gives second order convergence. For small values of ε ≪ h, the scheme is first order uniformly convergent in spatial direction.

Theorem 18.
Let u and U be the solution of (2)-(4) and (39), respectively. Then, the following uniform error bound holds: Proof. The combination of temporal and spatial error bounds gives the required result.

Numerical Examples, Results, and Discussions
Here, we illustrate the proposed scheme using model examples. The exact solutions of the considered examples are not known. We investigate the theoretical results by performing experiments using the proposed scheme.
Example 19. Consider singularly perturbed time delay parabolic PDEs: Example 20. Consider singularly perturbed time delay parabolic PDEs: with interval condition uðx, tÞ = 0, on ðx, tÞ ∈ ½0, 1 × ½−1, 0 and the boundary conditions Since and the ε-uniform (parameter uniform) error estimate by using We calculate the rate of convergence of the proposed scheme using and the ε-uniform rate of convergence using

Advances in Mathematical Physics
In Figures 1 and 2, the numerical solution of Examples 19 and 20 is given, respectively, for different values of perturbation parameterε. As we observe on these figures, a strong boundary layer is maintained on the right side of the spatial domain as ε goes small starting from ε = 2 −5 to 2 −20 . In these figures, we observe that the proposed scheme approximates the solution without creating oscillations or divergence. Note that time delay does not have an effect on position and size of the boundary layer since the layer occurs on the spatial domain direction. If the delay was on the spatial variable, we may have an interior layer which is created because of the delay parameter. In Tables 1 and 2, the maximum absolute error of Examples 19 and 20 is given, respectively, for different values of perturbation parameter ε and mesh numbers N, M. As we observed the results in these tables for each number of mesh N, M as ε becomes small, the maximum absolute error after showing growth becomes stable and identical. This indicates the stability and uniform convergence of the scheme. In the last two rows of these tables, the ε-uniform error and ε-uniform rate of convergence of the scheme are given. As one observes, the scheme gives first order uniform convergence. The second order convergence of the temporal discretization is depicted in Table 3. In Tables 4 and 5, comparison of the ε-uniform error and ε -uniform rate of convergence of the proposed scheme with results of some recently published papers is given. As we observed, the proposed scheme gives a more accurate result than the results in [16,[18][19][20] and [24]. The proposed scheme has a limitation for solving nonlinear singularly perturbed problems.

Conclusion
In this paper, singularly perturbed parabolic convectiondiffusion-reaction equation with large time delay is considered. The solution of the considered problem exhibits boundary layer of thickness ε on the right side of the spatial domain. The bounds and properties of the analytical solution are discussed. To stabilize the influence of the singular perturbation on the discrete solution, and exponential fitting parameter is developed using the derivation from the asymptotic solution. The numerical scheme is developed using the Crank-Nicolson method in temporal discretization and an exponentially fitted central finite difference method on the spatial discretization. The existence of the discrete solution is discussed using the comparison principle. The stability and uniform convergence of the scheme are investigated well theoretically. Numerical results are depicted using maximum absolute error, ε-uniform error, and ε-uniform rate of convergence in tables which are in good agreement with the theoretical analysis. The developed scheme gives stable and uniformly convergent result with linear order of convergence.

Data Availability
No additional data is used in this manuscript.

Conflicts of Interest
The authors declare no conflict of interest.