Abstract

A numerical scheme is developed to solve a large time delay two-parameter singularly perturbed one-dimensional parabolic problem in a rectangular domain. Two small parameters multiply the convective and diffusive terms, which determine the width of the left and right lateral surface boundary layers. Uniform mesh and piece-wise uniform Shishkin mesh discretization are applied in time and spatial dimensions, respectively. The numerical scheme is formulated by using the Crank–Nicolson method on two consecutive time steps and the average central finite difference approximates in spatial derivatives. It is proved that the method is uniformly convergent, independent of the perturbation parameters, where the convection term is dominated by the diffusion term. The resulting scheme is almost second-order convergent in the spatial direction and second-order convergent in the temporal direction. Numerical experiments illustrate theoretical findings, and the method provides more accurate numerical solutions than prior literature.

1. Introduction

Most mathematical models of real-life phenomena in physics, chemistry, biology, astronomy, and other fields of applied science involve the dependent variable and its derivatives as arguments in the equation. In many cases, the higher-order terms are multiplied by “small” interrelated multiparameters, which are referred to as singularly perturbed differential equations. If we consider two perturbation parameters, one (say μ) can be assumed to be a function of the other (say ). The earliest and popular paper [1] introduces these problems as two-parameter singularly perturbed differential equations defined aswhere and are “small” positive interrelated parameters simultaneously approaching zero and , , and are linear differential operators whose orders are , respectively, for a sufficiently smooth function . These kinds of equations arise in the fields of fluid dynamics, quantum mechanics, chemical flow reactor theory, and DC motor analysis [2, 3].

The two-parameter singularly perturbed parabolic differential equation is the focus of some recent studies [411]. Some others considered this problem with nonsmooth data, which implies boundary and interior layers in their solution [1216]. However, they did not consider the delay condition. In other models of physical problems like heat and mass transfer, control theory, and chemical processes, the system may depend on the present and past history of the state which yields a delay term in the differential equation. Such singularly perturbed time delay parabolic problems with uniperturbation parameter are addressed by some literature [1723]. Besides solving singularly perturbed differential equations, obtaining higher-order and robust numerical solutions is the main devotion of researchers [2430].

The concern of our problem is a two-parameter singularly perturbed time delay one-dimensional parabolic differential equation for the given initial and two boundary conditions. Different numerical methods of this type of problems are studied by some authors [3134]. Govindarao and others in [31] solved the problem by applying implicit Euler finite difference approximation on uniform mesh in the direction of time and upwind finite difference approximation on Shishkin and Bakhvalov–Shishkin mesh in space direction. They proved that their scheme is parameter-uniform with where and are spatial and temporal mesh sizes. Kumar and Others in [32] improved the spatial order of convergence to almost second by implementing a hybrid monotone finite difference scheme in the direction of space. Recently, two articles entitled “An effective numerical approach for two-parameter time-delayed singularly perturbed problems” and “Fitted cubic spline scheme for two-parameter singularly perturbed time-delay parabolic problems” also addressed this problem. The first article applied the Crank–Nicolson scheme in temporal direction and the B-spline collocation method for spatial [33] and the second article applied a combination of the fitted operator scheme of a cubic spline with —method on a uniform mesh grid [34].

The properties of the solution to these problems mainly depend on the dominant of diffusive or convective parameter. The problem has a reaction-diffusion problem property and is referred to as diffusion-dominant if the convective parameter is significantly smaller than the diffusive parameter. This paper proposes to develop a more accurate and parameter-uniform numerical scheme to solve a two-parameter singularly perturbed one-dimension delay parabolic problem using the Crank–Nicolson finite difference approximation on a uniform mesh in the time direction and the central finite difference method on Shishkin mesh in the spatial direction when the diffusive parameter dominates the convective parameter.

This article consists of the following structures: in Section 2, the governing problem will be stated and the appropriate regularity of the analytical solution will be discussed. Section 3 discusses some prior bounds of the analytical solution and its derivatives. The discretization of the domain and development of the numerical scheme are presented in Section 4. Section 5 contains the stability and convergence analysis of the scheme. Section 6 strengthens theoretical results using numerical solutions of counterexamples with tables and graphs, and finally, a brief conclusion is given in Section 7.

In every part of this article, denotes the supreme norm, if is a continuous function defined on , and if is a discrete value on the tensor of the domain, where is the discrete parameter. denotes the generic positive constant independent of perturbation parameters and mesh sizes. The Landau symbol is used to indicate order relations.

2. Statement of the Problem

We consider the two-parameter singularly perturbed large time delay one-dimensional parabolic problem defined on the rectangular domaingiven bywith initial and boundary conditionswhere and are perturbation parameters such that and , the delay parameter , such that the terminal time, for some positive integer . The functions and are sufficiently smooth and bounded that satisfyassuming that and are sufficiently smooth and bounded on the boundary, and compatible at the corner points.

For , the problem is treated as a standard two-parameter singularly perturbed parabolic problem since the value is known. The computation extends for after solving for . Repeating this method of steps completes the solution for the whole time domain component.

The problem is a kind of reaction-diffusion when and convection-diffusion when , but our interest is in the case . The analysis was made concerning two cases of O’Malley’s study in [1], when and . In the first case, the problem resembles the properties of the case when , and hence, layers appear in the neighborhood of and . For the second case, layer widths and appear in the neighborhood of and , respectively.

3. Bounds of Continuous Solution and Its Derivatives

Let us define a differential operator on the continuous solution of problems (3) and (4) by

It is assumed that the data are compatible at its boundary and then the following compatibility conditions should satisfy at the corner points , , , and .which makes sense whenever , , and . The uniqueness and existence of are established depending on these conditions, and it is also necessary to show the differential operator, and satisfies the maximum principle.

Lemma 1. Continuous maximum principle
Suppose , such that and . Then, .

Proof. Refer [34]

Lemma 2. The solution of equations (3) and (4) satisfies the following bounds:for some constant independent of and .

Proof. Set two barrier functions defined asapplying the differential operatorFor sufficiently large C, we haveSince satisfies the maximum principle, it holds the comparison relation .
Using the relation , as T is the maximum of t, , then .
Taking the supreme of all , we arrive .

Lemma 3. For non-negative integers i and j, such that , the derivatives of solution of problems (3) and (4) satisfy the following bounds:

Proof. (see [6]).
For the purpose of error analysis, it is necessary to decompose the solution of problems (3) and (4) into a regular component , which characterizes the solution behavior outside the boundary layers, and a singular component , which characterizes the solution behavior inside the boundary layers. Further, we decompose the singular component into the sum of left and right singular parts such that . The components individually satisfySome prior bounds of the components of the continuous solution and their derivatives are given in the following lemmas, which are necessary for the later error estimate.

Lemma 4. The bounds of the left and right singular components have given by

Proof. (see [6]).

Lemma 5. The derivatives of regular and singular components of satisfy the following bounds:(i)If then(ii)If then

The derivative bounds of and satisfy the derivative bounds of u (x, t) under Lemma 3.

Proof. Refer [32].

4. Formulation of Numerical Scheme

Applying the Crank–Nicolson approach upon uniform discretization for time and using midpoint central finite difference approximation on Shishkin mesh type in spatial direction, the derivation is explained in the following subsections.

4.1. Semidiscretization

The temporal discretization is performed by uniform mesh extending up to the time-lag interval aswhere M and m are the number of meshes in the intervals and , respectively, assuming that M is a positive multiple of m.

For problems (3) and (4) at , the mean point of and , for all , is given by

Applying Taylor’s series expansion of and about we have

From these expansions, we can derive

At time step, we substitute (22) for time derivative and average values between and for spatial derivatives into (20) and we obtainwhere is the identity differential operator, , and that possess a truncation error of .

Equation (23) gives the semidiscrete boundary value problem formulation of problems (3) and (4) for each time step with the boundary conditions as

Let be the numerical solution of the semidiscrete problem whose exact solution is at a fixed and is the truncation error caused by one step of iteration which is called local truncation error. Hence, the local absolute maximum error . The global truncation error is the cumulative truncation error up to iteration. Let E determine (the global error throughout the whole time interval).

Lemma 6. The global error estimate, of the semidiscrete problem (23) and (24), is given by

Proof. The local maximum absolute error, , is bounded byThus,
Considering the semidiscretization, at each time step , the differential operator , defined in 23 satisfies the maximum principle.

Lemma 7. Semidiscrete maximum principle
For a fixed time , let the function , satisfies and . Then, .

Proof. We apply a similar technique of contradiction as the continuous maximum principle of Lemma 1.
Let as . Supposition of implies a contradiction of the given hypothesis and hence proved.

Lemma 8. Semidiscretization uniform stability
The solution of equations (23) and (24) satisfies

Proof. Define barrier functionsWe haveThus, applying Lemma 7 gives the required result.
The boundary layer properties are determined according to the characteristic equation of (23), which iswhose roots areLet

Lemma 9. The bound of the solution, of the discrete problem (23) and (24), is provided byup to a certain prescribed order q, for and .

Proof. Refer [5, 35].

4.2. Full Discretization

Depending on the prior knowledge of boundary layers’ width and position, we apply piece-wise uniform Shishkin mesh type for spatial discretization. Let the transition points be and . Then, divide the spatial domain component into three subintervals , and . According to Shishkin mesh [36], each subinterval and will be divided into a uniform mesh size of and , respectively, where N is the number of meshes in spatial direction assuming that it is a multiple of 4. The transition points and are determined by

So, the spatial discretization, where

Extending the semidiscretized equation (23) to full discretization using central difference approximation for a fixed time , since the spatial discretization is nonuniform mesh, we have the finite difference approximations,

and for each mesh point Note that . Then, the finite difference scheme of equations (23) and (24) becomesor equivalently, using the Crank–Nicolson finite difference formula,for all and

Defining a difference operator , the corresponding difference equation can be written aswhich is

Taking all terms with unknown nodal values to the left and known values to the right, it has a recurrence form of

For each time step , for all , we drive an tridiagonal system of equation,given all the discrete initial and boundary valueswhere

5. Convergence of Numerical Method

Lemma 10. If a square matrix is real, irreducible, diagonally dominant with , and for all and , then where is a zero matrix.

Proof. (see[37]).
This lemma can be a useful tool for establishing the conditions under which the scheme’s coefficient matrix has an inverse. Let us examine each hypothesis in turn.(i)The coefficient matrix of the system (41) is real,(ii)As , all diagonal entries are positive,(iii)Under the assumption of , for coarse mesh size , the lower and upper diagonal entries, and . Then, all off-diagonal entries of the coefficient matrix are nonpositive,(iv) and , which implies Since all the rest entries in row are zero, the matrix is diagonally dominant(v)The irreducibility of the matrix can be shown using directed graph analysis. If the directed graph of a matrix is strongly connected, then it is irreducible. One can refer [37] for more about the directed graph. For the positional nonzero entry inhibiting the self-mapping of the diagonal entries, the directed graph of the tridiagonal coefficient matrix resembles Figure 1, which is strongly connected, and hence, the coefficient matrix is irreducible.The descriptions through satisfy all the hypotheses of Lemma 10, and then the coefficient matrix is invertible, and we guarantee to determine the numerical solution uniquely. These properties are also verification of the matrix A being an M-matrix, which is a monotone matrix. The monotonicity of a matrix establishes the discrete maximum principle.

Lemma 11. Discrete maximum principle
Given a discrete function, with and and . If and , then and

Remark 12. It is also necessary to investigate the assumption for the two major cases and . In both cases, is greater than or equal to . This comparison with the assumption implies which isConsidering the first case, , there exist such that the assumption (44) holds for all . More generally, when as , the system (41) converges for all for some positive integer . But for the other case , the assumption 34 holds for some fortunate large N. That is, when as , the system does not converge.

5.1. Error Bounds

Alike the continuous solution decomposed into regular and singular components in Section 3, , the numerical solution can be also decomposed as with a similar property, which is discretely defined asthen the error

The triangle inequality law gives its bound with

As we discuss in Section 2, the numerical solution is intended to compute in the time step of one after the other. So, the error analysis is performed in two intervals of time, when or , and or .Case I: When The source function, in the right-hand side, the retarded term as well, is a known function. For this case, the analysis is computed as a nondelay differential equation.

Lemma 13. The singular components of the discrete solution satisfy

Proof. (see [32]).
The evaluation of the error estimate of each component starts with the truncation error of the components separately. The continuous regular component satisfies the continuous problem (13), and the discrete regular component satisfies the discrete problem (45). From differential equation (6) and difference equation (38), we haveHence, the truncation bound of the regular component isSimilarly, the singular components satisfy the continuous problem (14) and (15) and discrete problem (46) and (47). Then, we deriveand

Lemma 14. The left layer component satisfies the following error bound estimate:

Proof. In the left boundary layer region, , , if and , if . Applying the derivative bounds, equation (??) of Lemma 5, the bound of equation (53) givesConsidering , when and , when , we arriveThe stability of discrete scheme impliesBut out of the left boundary layer region, , the absolute error can be calculated using the triangle inequality rule and the bounds of individual components in Lemmas 4 and 13.Then, taking the supreme, it impliesEquations (58) and (60) give the result.

Lemma 15. The right layer component satisfies the following error bound estimate:

Proof. In the right boundary layer region, , if , then . Otherwise, . The derivative bounds of Lemma 3, singular components bounds of Lemma 4, Lemma 13, and the bound of equation (54) with a similar argument to that of Lemma 14 implies the result.

Lemma 16. The regular component satisfies the following error bound estimate:

Proof. In the left boundary layer region, the uniform mesh size, if or if Applying the derivative bounds of Lemma 3, the bound of equation (60) becomesConsidering the right boundary layer region, the uniform mesh size, if or if We write the bound of equation (60) asThe outer boundary layer region, the mesh size, , thenCombining the three spatial interval cases and applying the maximum principle, we arrive at the required result.Case II: When In this case, the approximated values are involved on the right-hand side. Then, the error propagated in these time intervals must be shown, not magnified. We consider the left singular components asBy Lemma 14,Similar arguments can be given for the right layer and outer layer bounds, then Lemma 14 through 16 holds for the whole domain  × [0, T). We conclude by the following theorem.

Theorem 17. The fitted mesh scheme (38) is uniformly convergent under the condition . The exact solution, of problem (3) and (4) and approximate solution, of problem (41) and (42), satisfies the following error estimate.

6. Numerical Results and Discussion

Two examples are adopted from [31] for the numerical experiment. Numerical solutions, maximum error estimates, convergence rates, and comparison with prior literature are demonstrated for particular perturbation parameters and mesh grids in graphs and tables.

Example 1. subjected to the initial and boundary conditions

Example 2. subjected to the initial and boundary conditions

Graphs of numerical solutions applying the developed scheme (41) and (42) for are shown in Figures 2 and 3.

The exact solutions of the considered examples are not known. So, the point-wise error estimate and rate of convergence of the proposed scheme are computed using the double mesh principle. For each perturbation parameter and , we set numerical solutions on a grid mesh and on a doubled grid mesh . The error estimate of the numerical solution on the discrete parameters and , isand the corresponding rate of convergence, is given as

The maximum absolute error, and the corresponding rate of convergence, when both and mutually tend to zero are defined by

The error estimate and rate of convergence of numerical solution of the test problems using the given scheme (41), with the considered conditions, fixing one of the perturbation parameters and reducing the other are tabulated in Tables 1 and 2. Table 3 illustrates the maximum error estimate and the corresponding rate of convergence of Example 1. A similar illustration is given through Tables 46 for Example 2. The maximum point-wise error decreases uniformly as the number of meshes increase irrespective of perturbation parameters’ values with the second order of convergence which is expected from our theoretical result in the previous section. The logarithmic scaled graph in Figure 4 strengthens the predicted convergence rate.

The numerical solution in [34] compared the error and convergence rate of problem 1, when and for mesh grid , with the prior literature results of [32, 33]. We also include our numerical out comes in comparison with these prior results in Table 7.

7. Conclusion

In this study, a numerical method is developed to resolve a one-dimensional parabolic problem with a long time delay and two singularly perturbed parameters. Diffusion and convection terms are multiplied by the perturbation parameters. But we take into account areas where the diffusion term predominates. In the time and spatial dimensions, uniform mesh and piece-wise uniform Shishkin mesh discretization are used, respectively. The Crank–Nicolson method is used to formulate the numerical scheme on two successive time steps, and average central finite differences are used to approximate the spatial derivatives. Under the premise of diffusion-dominant problems , we computationally showed that the approach is uniformly convergent, independent of perturbation parameters. We also validated the numerical solutions of the governing equation by solving two test problems. The maximum point-wise error decreases uniformly independent of perturbation parameters with almost the second order of convergence as the number of meshes doubles, which aligns with the theoretical analysis. The logarithmic scaled graphs in Figure 4 demonstrate the same convergence rate. Although we restrict the method to convective-dominated problems, the formulated scheme is more accurate and efficient compared to numerical methods that exist in prior literature. We focus on diffusion-dominant and one-dimensional problems; one can extend this work to two-dimensional problems and study the other case, conviction-dominant problems.

Data Availability

No data were used to support the findings of this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors want to thank Arba Minch University for the material and financial support.