Abstract
A fitted tension spline numerical scheme for a singularly perturbed parabolic problem (SPPP) with time delay is proposed. The presence of a small parameter as a multiple of the diffusion term leads to the suddenly changing behaviors of the solution in the boundary layer region. This results in a challenging duty to solve the problem analytically. Classical numerical methods cause spurious nonphysical oscillations unless an unacceptable number of mesh points is considered, which requires a large computational cost. To overcome this drawback, a numerical method comprising the backward Euler scheme in the time direction and the fitted spline scheme in the space direction on uniform meshes is proposed. To establish the stability and uniform convergence of the proposed method, an extensive amount of analysis is carried out. Three numerical examples are considered to validate the efficiency and applicability of the proposed scheme. It is proved that the proposed scheme is uniformly convergent of order one in both space and time. Further, the boundary layer behaviors of the solutions are given graphically.
Keywords: delay differential equation; fitted tension spline scheme; singularly perturbed; uniform convergence
1. Introduction
Delay differential equations in which its higher order derivative term is multiplied by a small perturbation parameter () and involving a delay term is known as singularly perturbed delay differential equations (SPDDEs). Physical phenomena such as gestation times, incubation periods, transport delays, and more can be represented by delays/lags [1]. SPDDEs have numerous applications in science and engineering, especially in biological, chemical, electronic, and transportation systems. Some of these include simulation of oil extraction from underground reservoirs [2], chemical processes, fluid flows, and water quality problems in river networks [3] and mechanical systems [4, 5]. A simplified example of a time-delayed mathematical model which is used in an automatic control system of a furnace to produce metal sheets is given as follows [6]: where denotes the temperature distribution in a metal sheet and and denote the heat source and the velocity with which the metal sheet is moving, respectively. Both and are dependent on the term . The time delay of fixed length is induced because of the finite speed of the controlling device. SPDDEs exhibit a boundary layer, where the solution changes rapidly as approaches zero. The boundary layer is a narrow region placed on the left or right side of the domain where the dependent variable undergoes very rapid [4]. Due to the existence of the boundary layer, solving SPDDEs using analytical or standard numerical schemes is a challenging task [7]. On the contrary, standard numerical methods cause spurious nonphysical oscillations in the numerical result unless an acceptable large number of mesh points are used; this is computationally not feasible [1, 8, 9]. To overcome these kinds of computational problems, one needs to look for a sound numerical method that converges uniformly regardless of .
Recently, a number of numerical methods have been conducted by different authors. Ayele, Tiruneh, and Derese [10] constructed a hybrid scheme of a proper combination of the midpoint upwind scheme in the outer layer region and the cubic spline method in the boundary layer region on the piecewise Shishkin mesh in the spatial direction. Das and Natesan [1] constructed a numerical scheme using the hybrid method on piecewise uniform Shishkin mesh in the space direction and an implicit Euler scheme in the temporal direction. Also, Govindarao and Mohapatra [11] constructed a numerical scheme using the hybrid scheme on Shishkin-type meshes in the space direction and an implicit trapezoidal scheme in the temporal direction. Gowrisankar and Natesan [12] proposed a numerical scheme using the upwind finite difference scheme on a piecewise uniform mesh. Salama and Al-Amery [13] constructed a numerical scheme comprising the Crank-Nicolson method for discretizing the time variable and the operator compact implicit (OCI) method in the spatial direction. Podila and Kumar [14] constructed a stable finite difference method that works on a uniform mesh and an adaptive mesh. The upwind finite difference method on Shishkin mesh is constructed in [15]. An adaptive mesh refinement approach is employed using the concept of entropy function in [16, 17]. An exponentially fitted finite difference method is considered in [9, 18]. The nonstandard finite difference scheme is considered in [19, 20] following Micken’s type discretization for the space derivatives. The authors in [7, 21] constructed the numerical techniques that work for large or small time delays. For the spatiotemporal delays differential equation, Ejere et al. [8] constructed a robust numerical technique. For the differential-difference equations, the authors in [22] developed the cubic spline in tension.
To the best of the authors’ knowledge, the fitted tension spline scheme has not been developed for solving singularly perturbed parabolic problems (SPPPs) with a time delay. As stated in [23], constructing uniformly convergent numerical methods regardless of is an active research area. This motivated us to look for a sound numerical method that converges uniformly regardless of . In this study, we developed a numerical scheme that comprises the backward Euler in the temporal direction and the fitted tension spline method in the space direction on uniform meshes. Thus, the main objective is to propose -uniform numerical method to treat SPPPs with a time delay.
In this study, the symbol denotes a positive constant that is independent of and the mesh elements. The norm , denoted by , is the maximum norm.
2. Continuous Problem
On the domain , we considered a SPPP of the form where .
and are given constants; the functions on and on are sufficiently smooth and bounded to satisfy the following assumptions: (i) . (ii) . Under assumption (i), the solution of Problem (2) exhibits a boundary layer of width along [1], and under assumption (ii), Problem (2) exhibits a boundary layer of width along [16].
2.1. A Priori Bounds
The existence and uniqueness of the solution of Equation (2) can be ensured by the sufficient smoothness of and , and the compatibility conditions of the corner points and the delay term [24] are stated below
Setting in Equation (2) obtained the reduced problem. Under assumption (i), the reduced problem has the form
Under assumption (ii), the reduced problem has the form
Lemma 1 (see [1]) (maximum principle). Let , given that and , then .
Lemma 2 (see [1, 11]). The solution of Equation (2) is estimated as
Lemma 3 (see [1, 11]). The solution of Equation (2) is estimated as
Lemma 4 (stability result). The solution of Equation (2) is estimated as where .
Proof 1. Determine the barrier functions as , on the boundaries, we get , , and Then,
which imply .
Thus, from Lemma 1, , such that,
Therefore, the proof is done.
Lemma 5 (see [1, 13]). Let be the solution of Equation (2), then its derivative satisfies the estimate for the boundary layer along and for the boundary layer along .
3. Numerical Method
3.1. The Time Semidiscretization
The time interval is divided uniformly with mesh points as . As stated in [25], the delay parameter is treated using the Taylor approximation. Using the backward Euler formula, we get
Equivalently, it is convenient to rewrite Equation (7) as where and with the boundaries
Equations (8) and (9) can be rewritten as with the boundaries where .
The local error estimate at the time step is given as
Lemma 6. If and , then is bounded as and the global error estimate
Proof 2. The proof is considered in [26].
Lemma 7. For , , then the derivative of the solution of Equations (8) and (9) satisfies the following bounds for the boundary layer along and for the boundary layer along .
Proof 3. For the proof, see [27].
3.2. The Spatial Discretization
The space interval divided into equally spaced nodes with mesh length is given as and .
3.2.1. Description of the Method
A function in that interpolates at the mesh points that depends on the compression parameter and is reduced to a cubic spline on as is referred to as a parametric cubic spline function [28]. For each interval , the spline function has the form where for is called cubic spline in tension. Putting and from the homogeneous part of Equation (14), we get for arbitrary constants and . Let the nonhomogeneous part be given as
Thus, we have where and . Using Equations (15) and (16), we obtain
From Equation (17), the arbitrary constants can be determined from the interpolation conditions and , then
The derivative of Equation (18) at on the interval is and on the interval is
Equating Equations (19) and (20) at the point , we get
Simplifying Equation (21), we obtain where and . The continuity condition in Equation (22) ensures the continuity of the first derivative of at interior nodes. Using into Equation (10), we obtain
The local truncation error obtained from Equation (22) is for any choice of and whose sum is 1/2. For the choice and , we have
Using Taylor series expansions, we approximate and as
Substituting Equation (24) into Equation (23) and substituting the resulting equation into Equation (22) and arranging, we get where .
3.2.2. Exponential Fitting Factor
Here, we find the exponential fitting factor to handle the effect of in the layer. As the theory of singularly perturbed problem given in [29] and Taylor’s series expansion of about , the zero-order asymptotic solution of Equations (10) and (11) is given as where .
Multiplying Equation (25) by with a term containing , we get
Thus, we consider two cases of the boundary layers:
Case 1. Left-end boundary layer for .
Multiplying Equation (27) by and taking a limit as , which gives
Substituting Equation (29) into Equation (28), we obtain
Case 2. Right-end boundary layer for .
Multiplying Equation (27) by and taking a limit as , which gives
Substituting Equation (32) into Equation (31), we obtain
Therefore, we take the variable exponential fitting factor
Therefore, inducing the obtained exponential fitting factor in Equation (34) into Equation (25) for and , we obtain where
In the explicit form, we write Equation (35) as where
3.3. Convergence Analysis
Lemma 8 (discrete comparison principle). There is a comparison function such that for and if and , then for .
Proof 4. For the proof, see [27, 30].
Lemma 9. The solution of Equation (35) satisfies where .
Proof 5. Let and the barrier functions defined as . On the boundaries, we obtain and . On the discretized space domain we have From Lemma 8, we obtain . Hence, the desired bound is obtained.
From Taylor’s series approximation, yield the bounds where .
For the constants and , we have
Following Equation (39), we have since .
The next theorem provides the error bound in the spatial direction for the boundary layer along .
Theorem 1. Let be the solution of Equations (8) and (9) and the solution of Equation (35) satisfies the following error estimate:
Proof 6. In the space direction, the error is given as Using the bounds in Equations (38) and (40), we obtain From Lemma 7, we obtain since and the desired bound is obtained.
Lemma 10. For a fixed mesh and as , it gives where and
Proof 7. The proof can be done by using L’Hospital’s rule. For the details, refer to [30].
Theorem 2. The solution of Equation (35) satisfies the following uniform error estimate:
Proof 8. From Lemma 10 and Equation (41), we get Hence, the result leads Using the over all , we get From Equation (42), we take two cases: when , the obtained method gives a second-order uniformly convergent. On the other hand, when , the method is first-order uniformly convergent in the space direction.
Theorem 3. Let and be the solutions of Equations (2) and (35), respectively. Then, the next uniform error estimate holds
Proof 9. The proof can be done following Lemma 6 and Theorem 2.
4. Numerical Results
The applicability of the proposed scheme is validated using three model examples. The numerical values are given for and . For the unknown exact solution of the model examples, we use the double mesh principle for the numerical experiment [8, 9]. Therefore, the maximum pointwise error , -uniform error , the rate of convergence , and uniform rate of convergence are calculated by , , and , respectively.
Example 1. Consider the problem [16] , with interval condition , on and the boundary conditions and , .
Example 2. Consider the problem [1] with interval condition on and the boundary conditions and ,
Example 3. Consider the problem [15] with interval condition on and the boundary conditions and ,
The , , and the corresponding of the proposed method are revealed in Tables 1–3 for each example, respectively, for different values of and . These tables show that for every value of , monotonically decreases as the step sizes decrease, and as , the after showing growth and remains constant, demonstrating -uniform convergence of the developed method. On the other hand, the evaluated and applying the developed method are displayed in the last two rows, which support the theoretical result of the proposed method in the first-order in the space direction.
The numerical solutions of the developed method for each example are revealed in Figures 1–3. From Figure 1, we verify that a strong boundary layer is maintained along as . Moreover, Figures 2 and 3 have revealed that a strong boundary layer is maintained along as . Further, to show the effect of on the steepness of the boundary layer, the solutions are depicted in Figure 4. In each figure, we observe that as , the width of the boundary layer decreases, which confirms the desired result, that is, the boundary layer width of . In addition, in Figure 5, the of the method is plotted by the log–log scale. From these figures, one can see that the decreases as the step size decreases for every value of , which confirms the -uniform convergence of the developed method. The comparison of the findings of the proposed scheme with the existing recently published works is revealed in Tables 4–6. As we have seen in Table 4, the results obtained in this paper have better accuracy than those of [18, 21]. In addition, in Tables 5 and 6, it is shown that the obtained results are more accurate compared to that of [21].
5. Conclusion
A fitted tension spline numerical scheme is constructed to solve a SPPP with a time delay. The proposed method considered the SPPP to exhibit a parabolic boundary layer in the neighborhood of the left or right side of the domain as approaches . The solution varies abruptly in the layer region due to the presence of the perturbation parameter. The problem behaves in a multiscale character where the solution has a rapid change in the boundary layer region and is uniform otherwise. Due to this, classical numerical methods on a uniform mesh lead to an oscillatory solution and thus fail to give satisfactory results. To overcome this drawback, we proposed a fitted spline numerical method. The method comprises the backward Euler in the temporal direction and the fitted spline method in the spatial direction on uniform meshes. The developed scheme gives an accuracy of first-order both in space and time direction. Three numerical examples are considered to validate the efficiency and applicability of the proposed scheme, and the results support our theoretical findings. Concisely, the presented scheme is stable, -uniform, and provides more accurate numerical results compared to those of others.
Data Availability Statement
All the generated data was used. No additional data were used to support the findings of this research work from other sources.
Conflicts of Interest
The authors declare no conflicts of interest.
Funding
The authors received no specific funding for this work.