A Numerical Approach for Diffusion-Dominant Two-Parameter Singularly Perturbed Delay Parabolic Differential Equations

. 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 difusive 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. Te numerical scheme is formulated by using the Crank–Nicolson method on two consecutive time steps and the average central fnite diference 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 difusion term. Te resulting scheme is almost second-order convergent in the spatial direction and second-order convergent in the temporal direction. Numerical experiments illustrate theoretical fndings, and the method provides more accurate numerical solutions than prior literature.


Introduction
Most mathematical models of real-life phenomena in physics, chemistry, biology, astronomy, and other felds 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 diferential equations.If we consider two perturbation parameters, one (say μ) can be assumed to be a function of the other (say ε).Te earliest and popular paper [1] introduces these problems as two-parameter singularly perturbed diferential equations defned as where ε and μ are "small" positive interrelated parameters simultaneously approaching zero and L 1 , L 2 , and L 3 are linear diferential operators whose orders are l 1 > l 2 > l 3 , respectively, for a sufciently smooth function φ.Tese kinds of equations arise in the felds of fuid dynamics, quantum mechanics, chemical fow reactor theory, and DC motor analysis [2,3].Te two-parameter singularly perturbed parabolic differential equation is the focus of some recent studies [4][5][6][7][8][9][10][11].Some others considered this problem with nonsmooth data, which implies boundary and interior layers in their solution [12][13][14][15][16].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 diferential equation.Such singularly perturbed time delay parabolic problems with uniperturbation parameter are addressed by some literature [17][18][19][20][21][22][23].Besides solving singularly perturbed diferential equations, obtaining higher-order and robust numerical solutions is the main devotion of researchers [24][25][26][27][28][29][30].
Te concern of our problem is a two-parameter singularly perturbed time delay one-dimensional parabolic diferential equation for the given initial and two boundary conditions.Diferent numerical methods of this type of problems are studied by some authors [31][32][33][34].Govindarao and others in [31] solved the problem by applying implicit Euler fnite diference approximation on uniform mesh in the direction of time and upwind fnite diference approximation on Shishkin and Bakhvalov-Shishkin mesh in space direction.Tey proved that their scheme is parameter-uniform with O(h + k) where h and k 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 fnite diference scheme in the direction of space.Recently, two articles entitled "An efective numerical approach for two-parameter timedelayed singularly perturbed problems" and "Fitted cubic spline scheme for two-parameter singularly perturbed timedelay parabolic problems" also addressed this problem.Te frst 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 ftted operator scheme of a cubic spline with θ-method on a uniform mesh grid [34].
Te properties of the solution to these problems mainly depend on the dominant of difusive or convective parameter.Te problem has a reaction-difusion problem property and is referred to as difusion-dominant if the convective parameter is signifcantly smaller than the diffusive parameter.Tis 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 fnite diference approximation on a uniform mesh in the time direction and the central fnite diference method on Shishkin mesh in the spatial direction when the difusive parameter dominates the convective parameter.
Tis 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.Te 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 fnally, a brief conclusion is given in Section 7.
In every part of this article, ‖.‖ denotes the supreme norm, ) is a continuous function defned on D, and ‖f i,j ‖ � max 1≤i≤N1≤j≤M |f(x i , t j )| if f i,j is a discrete value on the tensor of the domain, where N × M is the discrete parameter.C denotes the generic positive constant independent of perturbation parameters and mesh sizes.Te Landau symbol O is used to indicate order relations.

Statement of the Problem
We consider the two-parameter singularly perturbed large time delay one-dimensional parabolic problem defned on the rectangular domain given by with initial and boundary conditions where ε and μ are perturbation parameters such that 0 < ε ≪ 1 and 0 < μ ≪ 1, the delay parameter τ > ε, μ, such that the terminal time, T � kτ for some positive integer k.Te functions a(x, t), b(x, t), c(x, t), and f(x, t) are sufciently smooth and bounded that satisfy assuming that c b (x, t), c l (t), and c r (t) are sufciently smooth and bounded on the boundary, Γ � Γ b ∪ Γ l ∪ Γ r and compatible at the corner points.
the problem is treated as a standard two-parameter singularly perturbed parabolic problem since the value u(x, t − τ) is known.Te computation extends for (x, t) ∈ [0, 1] × [τ, 2τ] after solving u(x, t) for [0, 1] × [0, τ].Repeating this method of steps completes the solution for the whole time domain component.
Te problem is a kind of reaction-difusion when μ � 0 and convection-difusion when μ � 1, but our interest is in the case 0 < μ ≪ 1. Te analysis was made concerning two cases of O'Malley's study in [1], when μ 2 ≤ Cε and μ 2 ≥ Cε.In the frst case, the problem resembles the properties of the case when μ � 0, and hence, O( � ε √ ) layers appear in the neighborhood of x � 0 and x � 1.For the second case, layer widths O(ε/μ) and O(μ) appear in the neighborhood of x � 0 and x � 1, respectively.

Bounds of Continuous Solution and Its Derivatives
Let us defne a diferential operator L ε,μ on the continuous solution u(x, t) 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 (0, 0), (1, 0), (0, − τ), and (1, − τ).

2
International Journal of Mathematics and Mathematical Sciences which makes sense whenever Te uniqueness and existence of u(x, t) ∈ C 2,1 (D) are established depending on these conditions, and it is also necessary to show the diferential operator, and L ε,μ satisfes the maximum principle.
Proof.Set two barrier functions defned as applying the diferential operator For sufciently large C, we have Since L ε,μ satisfes the maximum principle, it holds the comparison relation |s(x, t)| ≤ q(x, t).i.e, |u(x, t) Taking the supreme of all (x, t) ∈D, we arrive ‖u(x, t)‖ ≤ C. □ Lemma 3.For non-negative integers i and j, such that 0 ≤ i + 2j ≤ 4, 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 u(x, t) of problems (3) and (4) into a regular component v(x, t), which characterizes the solution behavior outside the boundary layers, and a singular component w(x, t), which characterizes the solution behavior inside the boundary layers.Further, we decompose the singular component into the sum of left and right singular parts w(x, t) � w L (x, t) + w R (x, t) such that u(x, t) � v(x, t) + w L (x, t) + w R (x, t).Te components individually satisfy International Journal of Mathematics and Mathematical Sciences Some 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. Te bounds of the left and right singular components have given by Proof (see [6]).

□
Lemma 5. Te derivatives of regular and singular components of u(x, t) satisfy the following bounds: Te derivative bounds of w L and w R satisfy the derivative bounds of u (x, t) under Lemma 3.

Formulation of Numerical Scheme
Applying the Crank-Nicolson approach upon uniform discretization for time and using midpoint central fnite diference approximation on Shishkin mesh type in spatial direction, the derivation is explained in the following subsections.

Semidiscretization.
Te temporal discretization is performed by uniform mesh extending up to the time-lag interval as where M and m are the number of meshes in the intervals [0, T] and [− τ, 0], respectively, assuming that M is a positive multiple of m.
For problems ( 3) and ( 4) at (x, t j+(1/2) ), the mean point of (x, t j ) and (x, t j+1 ), for all x ∈ D x , is given by Applying Taylor's series expansion of (x, t j ) and (x, t j+1 ) about (x, t j+(1/2) ), we have From these expansions, we can derive At (p + 1) st time step, we substitute (22) for time derivative and average values between (x, t p ) and (x, t P+1 ) for spatial derivatives into (20) and we obtain 4 International Journal of Mathematics and Mathematical Sciences where I is the identity diferential operator, Equation (23) gives the semidiscrete boundary value problem formulation of problems ( 3) and ( 4) for each time step t p+1 with the boundary conditions as Let u j (x) be the numerical solution of the semidiscrete problem whose exact solution is u(x, t j ) at a fxed t � t j and is the truncation error caused by one step of iteration which is called local truncation error.Hence, the j th local absolute maximum error 3 .Te global truncation error TE j is the cumulative truncation error up to j th iteration.Let E determine E M (the global error throughout the whole time interval).Lemma 6. Te global error estimate, E of the semidiscrete problem ( 23) and (24), is given by Proof.Te j th local maximum absolute error, e j , is bounded by Tus, ‖E‖ ≤ C(Δt) 2 .

□ Lemma 8. Semidiscretization uniform stability
Te solution u(x, t p+1 ) of equations ( 23) and ( 24) satisfes We have International Journal of Mathematics and Mathematical Sciences Tus, applying Lemma 7 gives the required result.Te boundary layer properties are determined according to the characteristic equation of (23), which is whose roots are Let □ Lemma 9. Te bound of the solution, u(x, t p+1 ) of the discrete problem ( 23) and ( 24), is provided by up to a certain prescribed order q, for 0 ≤ i ≤ q and 0 < ρ < 1.
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 σ 1 and σ 2 .Ten, divide the spatial and and [1 − σ 2 , 1] will be divided into a uniform mesh size of (N/4), (N/2) and (N/4), respectively, where N is the number of meshes in spatial direction assuming that it is a multiple of 4. Te transition points σ 1 and σ 2 are determined by So, the spatial discretization, Extending the semidiscretized equation (23) to full discretization using central diference approximation for a fxed time t � t p , since the spatial discretization is nonuniform mesh, we have the fnite diference approximations, for each mesh point (x i , t p ).Note that u i,j � u(x i , t j ).Ten, the fnite diference scheme of equations ( 23) and (24) becomes 6 International Journal of Mathematics and Mathematical Sciences or equivalently, using the Crank-Nicolson fnite diference formula, for all i � 1, 2, . . ., N + 1 and j � 1, 2, . . ., M + 1. Defning a diference operator , the corresponding diference equation can be written as which 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 (j + 1), for all i � 2, • • • , N, we drive an (N − 1) × (N − 1) tridiagonal system of equation, given all the discrete initial and boundary values where International Journal of Mathematics and Mathematical Sciences

Convergence of Numerical Method
Lemma 10.If a square matrix A � [a i,j ] is real, irreducible, diagonally dominant with a i,j ≤ 0, and a i,i > 0 for all i and j, then A − 1 > 0, where 0 is a zero matrix.
Proof (see [37]).Tis lemma can be a useful tool for establishing the conditions under which the scheme's coefcient matrix has an inverse.Let us examine each hypothesis in turn.
(i) Te coefcient matrix of the system (41) is real, (ii) As > 0, all diagonal entries are positive, (iii) Under the assumption of ‖a‖μh < 2ε, for coarse mesh size h, the lower and upper diagonal entries, Ten, all of-diagonal entries of the coefcient matrix are nonpositive, (iv) Since all the rest entries in i th row are zero, the matrix is diagonally dominant (v) Te 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 i th nonzero entry P i , inhibiting the self-mapping of the diagonal entries, the directed graph of the tridiagonal coefcient matrix resembles Figure 1, which is strongly connected, and hence, the coefcient matrix is irreducible.
Te descriptions i through v satisfy all the hypotheses of Lemma 10, and then the coefcient matrix is invertible, and we guarantee to determine the numerical solution uniquely.Tese properties are also verifcation of the matrix A being an M-matrix, which is a monotone matrix.Te monotonicity of a matrix establishes the discrete maximum principle.
Remark 12.It is also necessary to investigate the assumption ‖a‖μh < 2ε for the two major cases αμ 2 ≤ ςε and αμ 2 ≥ ςε.In both cases, h is greater than or equal to 1/N.Tis comparison with the assumption implies (1/N) ≤ (2ϵ/αμ) which is Considering the frst case, αμ 2 ≤ ςε, there exist N o such that the assumption (44) holds for all N ≥ N o .More generally, when (μ 2 /ε) ⟶ 0 as ε ⟶ 0, the system (41) converges for all N ≥ N o for some positive integer N o .But for the other case αμ 2 ≥ ςε, the assumption 34 holds for some fortunate large N. Tat is, when (ε/μ 2 ) ⟶ 0 as μ ⟶ 0, the system does not converge.

Error Bounds.
Alike the continuous solution decomposed into regular and singular components in Section 3, (u � v + w L + w R ), the numerical solution can be also decomposed as U � V + W L + W R , with a similar property, which is discretely defned as International Journal of Mathematics and Mathematical Sciences then the error Te 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 0 < j < m or t ∈ [0, τ], and m + 1 Case I: When t ∈ [0, τ] Te source function, (f(x, t)) in the right-hand side, the retarded term as well, is a known function.For this case, the analysis is computed as a nondelay diferential equation.

Lemma 13. Te singular components of the discrete solution satisfy
(50) Proof.(see [32]).Te evaluation of the error estimate of each component starts with the truncation error of the components separately.Te continuous regular component satisfes the continuous problem (13), and the discrete regular component satisfes the discrete problem (45).From diferential equation ( 6) and diference equation (38), we have (51) Hence, the truncation bound of the regular component is Similarly, the singular components satisfy the continuous problem ( 14) and ( 15) and discrete problem ( 46) and (47).Ten, we derive and Proof.In the left boundary layer region, [0, Applying the derivative bounds, equation (??) of Lemma 5, the bound of equation ( 53) gives Te stability of discrete scheme implies But 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.
Ten, taking the supreme, it implies Equations ( 58) and (60) give the result.

□
Lemma 15.Te right layer component satisfes the following error bound estimate: Proof.In the right boundary layer region, Te 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. Te regular component satisfes the following error bound estimate:
International Journal of Mathematics and Mathematical Sciences Proof.In the left boundary layer region, [0, σ 1 ], the uniform mesh size, Applying the derivative bounds of Lemma 3, the bound of equation ( 60) becomes Considering the right boundary layer region, We write the bound of equation (60) as Te outer boundary layer region, International Journal of Mathematics and Mathematical Sciences Combining the three spatial interval cases and applying the maximum principle, we arrive at the required result.

Case II: When t ∈ (τ, T)
In this case, the approximated values are involved on the right-hand side.Ten, the error propagated in these time intervals must be shown, not magnifed.We consider the left singular components as By 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, 1] × [0, T).We conclude by the following theorem.

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.
Te 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  International Journal of Mathematics and Mathematical Sciences double mesh principle.For each perturbation parameter ε and μ, we set numerical solutions U i.j on a grid mesh M × N and U i,j on a doubled grid mesh 2M × 2N.Te error estimate of the numerical solution on the discrete parameters N and M, E N,M ε.μ is and the corresponding rate of convergence, R M,N ε.μ is given as  Te maximum absolute error, E N,M , and the corresponding rate of convergence, R M,N , when both ε and μ mutually tend to zero are defned by   Te error estimate and rate of convergence of numerical solution of the test problems using the given scheme (41), with the considered conditions, fxing 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 4-6 for Example 2. Te maximum pointwise 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.Te logarithmic scaled graph in Figure 4 strengthens the predicted convergence rate.
Te numerical solution in [34] compared the error and convergence rate of problem 1, when ε � 10 − 4 and μ � 10 − 9 for mesh grid 64 ≤ N � 4M ≤ 512, with the prior literature results of [32,33].We also include our numerical out comes in comparison with these prior results in Table 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.Difusion and convection terms are multiplied by the perturbation parameters.But we take into account areas where the difusion term predominates.In the time and spatial dimensions, uniform mesh and piece-wise uniform Shishkin mesh discretization are used, respectively.Te Crank-Nicolson method is used to formulate the numerical scheme on two successive time steps, and average central fnite diferences are used to approximate the spatial derivatives.Under the premise of difusiondominant problems (αμ 2 ≤ ςε), 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.Te 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.Te 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 efcient compared to numerical methods that exist in prior literature.We focus on difusion-dominant and one-dimensional problems; one can extend this work to two-dimensional problems and study the other case, conviction-dominant problems.
) Figure 1: Te directed graph of the coefcient matrix.
Bold values emphasize that it concludes the rate of convergence for any small parameter.

Table 6 :
Te maximum error estimate (E N,M ) and rate of convergence (R N,M ) of Example 2, for 10 − 40 ≤ ε, μ ≤ 10 − 16 .Bold values emphasize that it concludes the rate of convergence for any small parameter.