Fitted computational method for singularly perturbed convection-diffusion equation with time delay

Abstract

A uniformly convergent numerical scheme is proposed to solve a singularly perturbed convection-diffusion problem with a large time delay. The diffusion term of the problem is multiplied by a perturbation parameter, ε. For a small ε, the problem exhibits a boundary layer, which makes it challenging to solve it analytically or using standard numerical methods. As a result, the backward Euler scheme is applied in the temporal direction. Non-symmetric finite difference schemes are applied for approximating the first-order derivative terms, and a higher-order finite difference method is applied for approximating the second-order derivative term. Furthermore, an exponential fitting factor is computed and induced in the difference scheme to handle the effect of the small parameter. Using the discrete maximum principle, the stability of the scheme is examined and analyzed. The developed scheme is parameter-uniform with a linear order of convergence in both space and time. To examine the accuracy of the method, two model examples are considered. Further, the boundary layer behavior of the solutions is given graphically.

1. Introduction

Delay differential equations (DDEs) are differential equations in which the evolution of the system is influenced by its past history. DDEs are called retarded types if the delay argument does not appear in the highest-order derivative term; otherwise, they are neutral types. DDEs play an important role in a variety of fields, including robotics, biosciences [1], economics, epidemiology and mechanics [2], fluid dynamics, reaction-diffusion equations [3], and population dynamics [4].

A singularly perturbed delay differential equation (SPDDE) is a delay differential equation in which its higher-order derivative term is multiplied by a small perturbation parameter (0 < ε ≪ 1) and contains at least one delay parameter on the term different from the highest derivative. In contrast to the magnitude of the delay parameter with the perturbation parameter, the delay is classified as a large delay or a small delay. If the magnitude of the delay parameter of the SPDDE is smaller than the perturbation parameter, then the equation is said to be a singularly perturbed delay differential equation with a small delay, whereas when the magnitude of the delay parameter is higher than the perturbation parameter, it is said to be a singularly perturbed delay differential equation with a large delay [5]. A singularly perturbed problem, which arises as a time delay, occurs in many application areas of science and engineering, for instance, in the simulation of oil extraction from underground reservoirs, chemical processes, fluid flows, water quality problems in river networks, and mechanical systems [6].

The presence of ε as a multiple of the higher-order derivative term causes a boundary layer. The boundary layer is an asymptotically narrow region located in the neighborhood of the endpoints of the domain, where the solution has a steep gradient as ε tends to zero [6]. With the rapidly changing behavior of the solution in the boundary layer, one encounters computational difficulties in treating a singularly perturbed problem using analytically or classical numerical schemes. On the contrary, classical numerical schemes lead to spurious non-physical oscillations in the numerical solution, unless an unacceptably large number of mesh points are considered, which leads to a massive computational cost [7]. In response to this, different authors have to look for sounding numerical schemes which converge uniformly regardless of ε.

Recently, the authors in [8], proposed the implicit-Euler scheme in the time direction and the central difference scheme in the space direction. The authors in [7, 9], proposed the implicit-Euler scheme in the time direction and the hybrid scheme by a proper combination of the midpoint upwind in the outer region and the central difference scheme in the inner region in the spatial direction on the Shishkin mesh. Moreover, this method is addressed in [10], for the two-parameter problem. In [11], the authors proposed the implicit-Euler scheme in the time direction and the hybrid scheme on a generalized Shishkin mesh in the spatial direction. Gowrisankar and Natisan in [12] developed the backward Euler scheme in time direction and the upwind finite difference scheme in the spatial direction using a piecewise uniform mesh. The implicit Euler scheme in the time direction and the upwind scheme in the spatial direction are considered in [13]. In [14], the implicit trapezoidal scheme in the time direction and the hybrid scheme by proper combination of the midpoint upwind in the outer region and the central difference in the inner region in the spatial direction are used.

The implicit Euler scheme in the time direction and the central difference scheme in the space direction are used in [4]. The extended cubic B-spline is considered in [15]. A domain decomposition method is considered in [16, 17]. The authors in [18, 19] proposed hybrid scheme on both Shishkin and Bakhvalov meshes. Podila and Kumar [20] proposed a new stable finite difference scheme on a uniform mesh and also on an adaptive mesh. The backward Euler scheme in the time direction and exponentially fitted difference method is considered in [21]. The Crank-Nicolson method in the time direction and a novel fitted finite difference scheme in spatial direction are proposed in [22]. The Crank-Nicolson method in the time direction and an exponentially fitted spline in the spatial direction are discussed in [23]. The implicit Euler scheme in the time direction and the non-standard finite difference method in the space direction are considered in [24]. In [25], the authors proposed Crank-Nicolson method in the time direction and the operator compact implicit (OCI) method on the Shishkin mesh in the space direction. The backward Euler in the time direction and method of line following Micken's type discretization for the space derivatives are used in [26]. Sahoo and Gupta [27] used higher-order difference with an identity expansion (HODIE) on a piecewise uniform mesh. A similar technique was also used in [28] for a coupled system of singularly perturbed problems. The authors [29, 30] proposed the numerical schemes that work for both cases when the delay term is large or small.

The main aim of this work was to develop a ε-uniform numerical scheme for the class of singularly perturbed convection-diffusion problem with a large time delay. The method comprises the backward Euler scheme in the time direction and an exponentially fitted higher-order finite difference scheme in the spatial direction. Error bound and uniform convergence of the developed scheme is investigated and proved. The proposed scheme gives more accurate, stable, and uniformly convergent results.

In this study, C has been considered as a generic positive constant, which does not depend on Δs, Δt, and ε. The maximum norm is denoted by ‖.‖, which is defined by ‖γ‖ = max_s,t∈Ω|γ(s, t)|.

2. Continuous problem

Let Ω = Ω_s × Ω_t = (0, 1) × (0, 𝕋] for 𝕋 > 0, we consider SPDDE of the form where 𝔏_εz(s, t) = −εz_ss(s, t) + β(s, t)z_s(s, t) + α(s, t)z(s, t).

Here, ε ∈ (0, 1] and δ > 0 are the perturbation parameter and the delay parameter, respectively. We pretended that the functions β(s, t), α(s, t), κ(s, t), γ(s, t) on and ψ_b(s, t), ψ_l(t), ψ_r(t) on η = η_l ∪ η_r ∪ η_b are smooth enough and bounded which meet α(s, t) ≥ ϖ > 0, κ(s, t) ≥ φ > 0, β(s, t) ≥ μ > 0 on . These conditions assure that problem (1) has a boundary layer near s = 1 [7].

2.1. A priori bounds

Under the premises that the data are Hölder continuous and satisfy the following compatibility conditions at the corner points and the delay terms [31], we confirm the existence and uniqueness of the solution of (1) These assumptions and conditions are fulfilled. Then, the problem (1) admits a unique solution [31].

Setting ε = 0, the reduced problem of (1) is given as where z₀(s, t) is the solution of the reduced problem.

Lemma 2.1. Let z(s, t) be the solution of (1). Then, we have where C does not depends on ε.

Proof: The proof is considered in [7].

The operator in (1) satisfies the next lemma.

Lemma 2.2. (Maximum principle). Let ν(s, t) ∈ C²(Ω) ∪ C⁰(η) satisfies ν(s, t) ≥ 0 (s, t) ∈ η. If 𝔏ν(s, t) ≥ 0, (s, t) ∈ Ω, then.

Proof: The proof is considered in [14].

Lemma 2.3. (Stability result). Let z(s, t) be the solution of (1). Then, we have where ϖ ≤ α(s, t).

Proof: The proof is considered in [22].

Lemma 2.4. The derivative of the solution z(s, t) of (1) with respect to s and t satisfywhere μ ≤ β(s, t).

Proof: The proof is considered in [7].

3. Numerical scheme

3.1. Temporal semi-discretization

The time domain [0, 𝕋] is discretized uniformly with step size Δt as and with M + 1 mesh points in [0, 𝕋] and j + 1 mesh points in [−δ, 0]. We have 𝕋 = rδ for some positive integer r.

Applying the backward Euler scheme for time derivative, we get Simplifying (8), we have where and with the boundaries Now, (9) rewrite as where and and .

The local truncation error in the time direction is given as

Lemma 3.1. The local error e_mat t_msatisfies the bound

Lemma 3.2. The global error E_mat t_msatisfies the boundProof: Using Lemma 3.1, the global error E_m bound at mth time step is given as

Lemma 3.3. For every m = 0(1)M − 1, the solution Z^m(s) of (9)-(10) satisfies the estimateProof: From (11), , where g = Q(s) − P(s, t_m)Z(s).

Now, we integrate twice and we obtain where

Using inequality and the bound Hence, |Z_p(s)| ≤ C. Here, . The boundary condition Z(1) = 0 yields C₁ = 0. Now, the constants C₁ and C₂ must satisfy Since B(s) is bounded on (0, 1), B(1) − B(y) ≤ C(1 − y). Then, It follows that . Hence, . Finally, implies that The proof is done for i = 1. For i > 1 follows by induction and repeated differentiation. For the details, refer [32].

3.2. Spatial discretization

We discretize the spatial domain [0, 1] into N equal number of sub-intervals with the length of h as 0 = s₀, s₁, ..., s_N = 1, and s_n = nh, n = 0(1)N. Consider a smooth function Z(s) in the interval [0, 1]. From Taylor's series approximation, we get Following a similar relation of (15), it holds From (16), we have Substituting (17) into (16) and simplifying, we obtain where

From (11), we draw where Using (19), we have From the Taylor series approximations of and , we get Substituting (21) into (20), we have From (18), we draw Substituting (22) into (23) and rearranging, we obtain where

3.2.1. Computing the exponential fitting factor

We introduce the exponential fitting factor σ to handle the effect of ε in the layer. From the singular perturbation theory stated in [33], the zero order asymptotic solution of the problem of the form is given as From Taylor's series, approximation for β(s) and α(s) restricting to their first terms about s = 1 is given as where Z₀(s) is the solution of reduced problem. Taking h → 0 and solving (28) at s_n−1, s_n, and s_n+1, we get where . Multiplying (24) by h and the term containing ε by σ and evaluating the limit of the resulting equation as h → 0, we get From (29), we have Substituting (31) into (30) and simplifying yields Then, we get the fitting factor σ Therefore, the required scheme is taken as where In the explicit form, it becomes where

3.3. Stability and uniform convergence analysis

Lemma 3.4. (Discrete maximum principle). Assume thatand, then.

Proof: Assume that there is k ∈ {0, 1, 2, …, N}, such that . Assume that and from the assumption, it is shown that k ∉ {0, 1}. So, we have and . Then, we get for k = 1(1)N − 1. So, the assumption is wrong. Therefore, and ∀n = 0(1)N.

Lemma 3.5. (Uniform stability result). Letbe the solution of (33), then we havewhere P(s_n, t_m) ≥ ζ > 0.

Proof: Let and define the barrier functions by . On the boundaries, we get ψ_l(t_m) ≥ 0 and . For s_n, n = 1(1)N − 1, we obtain By Lemma 3.4, we get . Therefore, the needed bound is obtained.

From Taylor's series expansion, we get where

The next theorem provides the truncation error estimate for the developed scheme.

Theorem 3.6. Let the coefficients α(s, t_m), β(s, t_m), and κ(s, t_m) of (9)-(10) be sufficiently smooth such that Z^m(s) ∈ C⁴[0, 1]. Then, the solutionof (33) satisfies the next boundProof: The error estimate in the spatial direction is given as where and .

For the constants C₁ and C₂, we have for ρ ≤ 1. For ρ → ∞, since which gives |ρ coth(ρ) − 1| ≤ C₁ρ.

Generally, ∀ρ > 0, we express as and we have From the bounds in (38), (40), and (42), we have By Lemma 3.3, we have Obviously, ε⁻³ ≥ ε⁻², then we draw thus, it gives the wanted bound.

Lemma 3.7. For a fixed mesh and as ε → 0, it holdswhere s_n = nh, 1 ≤ n ≤ N − 1.

Proof: The proof is in [22].

Theorem 3.8. Letbe the solution of (33), then we have the following uniform error boundProof: Substituting Lemma 3.7 into (39), we arrive at

Hence, the result leads Using the sup over all ε ∈ (0, 1], we get From (46), when ε > h, the obtained method uniformly converges uniformly with order two in the space direction. When ε ≪ h, the method converges uniformly with order one in the space direction.

Theorem 3.9. Let z and Z are the solutions of (1) and (33), respectively, then we have the following uniform error boundProof: The proof is considered by combining of Lemma 3.1 and Theorem 3.8.

4. Numerical results

Considering two test examples we carry out some numerical inquiries to confirm the developed scheme is ε-uniform convergent. Since the exact solution of the examples are not known, we used a variant of double mesh principle is applied for the numerical inquiries. So, we calculate the maximum pointwise error by , the ε-uniform error by , the rate of convergence by and the ε-uniform rate of convergence by r^N,M = log2(E^N,M/E^{2N, 2M}).

4.1. Example

Consider the problem [7]

(s, t) ∈ (0, 1) × (0, 2] with interval condition z(s, t) = 0, on (s, t) ∈ [0, 1] × [−1, 0] and the boundary conditions z(0, t) = 0 and z(1, t) = 0, t ∈ [0, 2].

4.2. Example

Consider the problem [13]

(s, t) ∈ (0, 1) × (0, 2] with interval condition z(s, t) = 0, on (s, t) ∈ [0, 1] × [−1, 0] and the boundary conditions z(0, t) = 0 and z(1, t) = 0, t ∈ [0, 2].

For distinguishable values of ε and N, the obtained results for the model Examples 4.1 and 4.2, respectively, , , and r^N,M of the developed scheme are delineated in Tables 1, 2. From these tables, one can observe that the maximum absolute error decreases as the step sizes decrease for every value of ε, and as ε approaches to zero, the maximum absolute error after getting large becomes constant, which displays ε-uniform convergence of the proposed scheme regardless of ε. On the other hand, the calculated E^N,M and the corresponding r^N,M using the proposed scheme are given in the last two rows, which confirms that the theoretical finding of the developed scheme is order one in both space and time direction.

In Figures 1, 2, the numerical solutions of the method for Examples 4.1 and 4.2 for different values of ε are given, respectively, for N = 80 and M = 40. Figure 3 displays the effect ε on the solutions profile of the developed scheme for Examples 4.1 and 4.2. From the figures, we see that a strong boundary layer is created on the right side of the spatial domain as ε close to zero. Furthermore, in Figure 4, the maximum point wise errors of the scheme is shown by the log-log scale plot. From these figures, one can observe that maximum absolute error decreases as the step sizes decrease for every values of ε, which confirm ε-uniform convergence of the proposed scheme.

Figure 1

Figure 2

Figure 3

Figure 4

In Table 3, the comparison with results of the developed method with the existing recently published studies of [23, 29] are given for Example 4.1. In Table 4, the comparison with results of the developed method with the existing number of recently published studies of [15, 24, 29, 30] are given for Example 4.2. As one follows, the developed scheme holds more accurate.

5. Conclusion

We have developed a numerical method for solving singularly perturbed parabolic convection-diffusion equation with a large time delay. The solution of the problem exhibits a boundary layer on the right side of the domain. The solution has a steep gradient in the layer region due to the presence of ε. In the rapidly changing behavior of the solution in the layer region, one encounters computational difficulties to find the solution using analytically or using classical numerical methods. To handle this effect, we developed method comprises of the backward Euler scheme in the time direction and an exponentially fitted higher order finite difference scheme in the spatial direction. Using comparison principle, the stability of the discrete scheme is analyzed. The stability and uniformly convergent of the method are discussed theoretically. Numerical results are delineated by applying maximum point wise error, ε-uniform error and ε-uniform rate of convergence in tables which are in acceptable agreement with the theoretical analysis. The developed method contributes more accurate, stable, and ε-uniform with a linear order of convergence in the spatial and in the time direction. The proposed scheme can be extended for singularly perturbed turning point problems.

References

• 1.

  BakerCTBocharovGARihanFA. A Report on the Use of Delay Differential Equations in Numerical Modelling in the Biosciences. Citeseer (1999).

  • Google Scholar

• 2.

  GopalsamyK. Stability and Oscillations in Delay Differential Equations of Population Dynamics. vol. 74. Springer Science and Business Media (1992).

  • Google Scholar

• 3.

  BestehornMGrigorievaEV. Formation and propagation of localized states in extended systems. Ann Phys. (2004) 516:423–31. 10.1002/andp.200451607-806

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4.

  GovindaraoLMohapatraJDasA. A fourth-order numerical scheme for singularly perturbed delay parabolic problem arising in population dynamics. J Appl Math Comput. (2020) 63:171–95. 10.1007/s12190-019-01313-7

  • CrossRef
  • Google Scholar

• 5.

  TianH. Numerical Treatment of Singularly Perturbed Delay Differential Equations. Manchester, UK: The University of Manchester (2000).

  • Pubmed Abstract
  • Google Scholar

• 6.

  AnsariARBakrSAShishkinGI. A parameter-robust finite difference method for singularly perturbed delay parabolic partial differential equations. J Comput and Appl Math. (2007) 205:552–66. 10.1016/j.cam.2006.05.032

  • CrossRef
  • Google Scholar

• 7.

  DasANatesanS. Uniformly convergent hybrid numerical scheme for singularly perturbed delay parabolic convection–diffusion problems on Shishkin mesh. Appl Math Comput. (2015) 271:168–86. 10.1016/j.amc.2015.08.137

  • CrossRef
  • Google Scholar

• 8.

  GowrisankarSNatesanS. A robust numerical scheme for singularly perturbed delay parabolic initial-boundary-value problems on equidistributed grids. Electron Trans Numer Anal. (2014) 41:376–95. 10.1016/j.cpc.2014.04.004

  • CrossRef
  • Google Scholar

• 9.

  GuptaVKumarMKumarS. Higher order numerical approximation for time dependent singularly perturbed differential-difference convection-diffusion equations. Num Methods Partial Differ Equ. (2018) 34:357–80. 10.1002/num.22203

  • CrossRef
  • Google Scholar

• 10.

  SumitKumarSKuldeepKumarM. A robust numerical method for a two-parameter singularly perturbed time delay parabolic problem. Comput Appl Math. (2020) 39:1–25. 10.1007/s40314-020-01236-1

  • CrossRef
  • Google Scholar

• 11.

  KumarSKumarM. High order parameter-uniform discretization for singularly perturbed parabolic partial differential equations with time delay. Comput Math Appl. (2014) 68:1355–67. 10.1016/j.camwa.2014.09.004

  • CrossRef
  • Google Scholar

• 12.

  GowrisankarSNatesanS. ε-Uniformly convergent numerical scheme for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2017) 94:902–21. 10.1080/00207160.2016.1154948

  • CrossRef
  • Google Scholar

• 13.

  DasANatesanS. Second-order uniformly convergent numerical method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2018) 95:490–510. 10.1080/00207160.2017.1290439

  • CrossRef
  • Google Scholar

• 14.

  GovindaraoLMohapatraJ. A second order numerical method for singularly perturbed delay parabolic partial differential equation. Eng Comput. (2018) 36:420–44. 10.1108/EC-08-2018-0337

  • CrossRef
  • Google Scholar

• 15.

  KumarDKumariP. A parameter-uniform scheme for singularly perturbed partial differential equations with a time lag. Num Methods Partial Differ Equ. (2020) 36:868–86. 10.1002/num.22455

  • CrossRef
  • Google Scholar

• 16.

  SinghJKumarSKumarM. A domain decomposition method for solving singularly perturbed parabolic reaction-diffusion problems with time delay. Num Methods Partial Differ Equ. (2018) 34:1849–66. 10.1002/num.22256

  • CrossRef
  • Google Scholar

• 17.

  RaoSCSKumarS. An almost fourth order uniformly convergent domain decomposition method for a coupled system of singularly perturbed reaction–diffusion equations. J Comput Appl Math. (2011) 235:3342–54. 10.1016/j.cam.2011.01.047

  • CrossRef
  • Google Scholar

• 18.

  KumarSKumarM. A second order uniformly convergent numerical scheme for parameterized singularly perturbed delay differential problems. Num Algor. (2017) 76:349–60. 10.1007/s11075-016-0258-9

  • CrossRef
  • Google Scholar

• 19.

  KumarSKumarM. Analysis of some numerical methods on layer adapted meshes for singularly perturbed quasilinear systems. Num Algor. (2016) 71:139–50. 10.1007/s11075-015-9989-2

  • CrossRef
  • Google Scholar

• 20.

  PodilaPCKumarK. A new stable finite difference scheme and its convergence for time-delayed singularly perturbed parabolic PDEs. Comput Appl Math. (2020) 39:1–16.

  • Google Scholar

• 21.

  TesfayeSKWoldaregayMMDinkaTGDuressaGF. Fitted computational method for solving singularly perturbed small time lag problem. BMC Res Notes. (2022) 15:318. 10.1186/s13104-022-06202-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22.

  WoldaregayMMAnileyWTDuressaGF. Novel numerical scheme for singularly perturbed time delay convection-diffusion equation. Adv Math Phys. (2021) 2021:6641236. 10.1155/2021/6641236

  • CrossRef
  • Google Scholar

• 23.

  NegeroNTDuressaGF. An efficient numerical approach for singularly perturbed parabolic convection-diffusion problems with large time-lag. J Math Model. (2022) 10:173–90. 10.22124/jmm.2021.19608.1682

  • CrossRef
  • Google Scholar

• 24.

  BabuGBansalK. A high order robust numerical scheme for singularly perturbed delay parabolic convection diffusion problems. J Appl Math Comput. (2021) 68:363–89. 10.1007/s12190-021-01512-1

  • CrossRef
  • Google Scholar

• 25.

  SalamaAAAl-AmeryDG. A higher order uniformly convergent method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2017) 94:2520–46. 10.1080/00207160.2017.1284317

  • CrossRef
  • Google Scholar

• 26.

  NegeroNTDuressaGF. A method of line with improved accuracy for singularly perturbed parabolic convection–diffusion problems with large temporal lag. Results Appl Math. (2021) 11:100174. 10.1016/j.rinam.2021.100174

  • CrossRef
  • Google Scholar

• 27.

  SahooSKGuptaV. Parameter robust higher-order finite difference method for convection-diffusion problem with time delay. Num Methods Partial Differ Equ. (in press). 10.1002/num.23039

  • CrossRef
  • Google Scholar

• 28.

  RaoSCSKumarS. Second order global uniformly convergent numerical method for a coupled system of singularly perturbed initial value problems. Appl Math Comput. (2012) 219:3740–53. 10.1016/j.amc.2012.09.075

  • CrossRef
  • Google Scholar

• 29.

  KumarDKumariP. A parameter-uniform numerical scheme for the parabolic singularly perturbed initial boundary value problems with large time delay. J Appl Math Comput. (2019) 59:179–206. 10.1007/s12190-018-1174-z

  • CrossRef
  • Google Scholar

• 30.

  NegeroNTDuressaGF. Uniform convergent solution of singularly perturbed parabolic differential equations with general temporal-lag. Iran J Sci Technol Trans A Sci. (2022) 46:507–24. 10.1007/s40995-021-01258-2

  • CrossRef
  • Google Scholar

• 31.

  ProtterMWeinbergerH. Maximum Principles in Differential Equations. Englewood Cliffs, NJ: Prentice-Hall Inc. (1967).

  • Google Scholar

• 32.

  KelloggRBTsanA. Analysis of some difference approximations for a singular perturbation problem without turning points. Math Comp. (1978) 32:1025–39.

  • Google Scholar

• 33.

  O'MalleyRE. Singular Perturbation Methods for Ordinary Differential Equations. vol. 89. New York, NY: Springer (1991).

  • Pubmed Abstract
  • Google Scholar