Parameter-uniformly convergent numerical scheme for singularly perturbed delay parabolic differential equation via extended B-spline collocation

This paper presents a parameter-uniform numerical method to solve the time dependent singularly perturbed delay parabolic convection-diffusion problems. The solution to these problems displays a parabolic boundary layer if the perturbation parameter approaches zero. The retarded argument of the delay term made to coincide with a mesh point and the resulting singularly perturbed delay parabolic convection-diffusion problem is approximated using the implicit Euler method in temporal direction and extended cubic B-spline collocation in spatial orientation by introducing artificial viscosity both on uniform mesh. The proposed method is shown to be parameter uniform convergent, unconditionally stable, and linear order of accuracy. Furthermore, the obtained numerical results agreed with the theoretical results.

1. Introduction

Singularly perturbed delay differential equations (SPDDEs) are differential equations that involve diffusion parameter, known as the perturbation parameter ε, and at least a delay term. SPDDEs arise in many practical phenomena, such as epidemiology [1], population growth [2], chemostat models [3], circadian rhythms [4], the respiratory system [5], and tumor growth [6]. The following is a typical example of SPDDEs, which is a mathematical model of the overall control system [7] which models a furnace used to process metal sheets.

$\frac{\partial u\left(x,t\right)}{\partial t}-\epsilon \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}=\upsilon \left(g\left(u\left(x,t-\tau \right)\right)\right)\frac{\partial u}{\partial x}+C\left[f\left(u\left(x,t-\tau \right)\right)-u\left(x,t\right)\right],$

defined on a one-dimensional spatial domain 0 < x < 1, where u(x, t) is the temperature distribution in metal sheet moving at a velocity υ depending on a prescribed spatial average of the time-delayed temperature distribution u(x, t−τ) and f represents a distributed temperature source function depending on u(x, t−τ). The spatial temperature distribution in the incoming and out coming material within the furnace is given as u(0, t) and u(1, t), respectively. A controlling device monitoring the current temperature dynamically adapts both υ and f. A fixed delay of length τ is induced by the finite speed of the controller. Another typical example of SPDDEs is the following logistic equation [8]:

$\frac{\partial u\left(x,t\right)}{\partial t}-\epsilon \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}=u\left(x,t\right)\left(1-u\left(x,t-\tau \right)\right),$

which arise in mathematical ecology for the evolution of a population with density u(x, t). A population u(x, t) depends on the population at an earlier time, t − τ rather than t. The delay τcan arise from a great variety of causes, such as duration of gestation, hatching period, and slow replacement of food supplies. Thus, u(x, t) depends on average past population u(x, t−τ). The initial and final spatial spread of the favored population are given by u(0, t) and u(1, t), respectively. A set of examples is available in Wu et al. [9] to illustrate the wide range of existing delay differential equation models.

When the highest order term coefficient ε tends to zero, the solution of SPDDE possesses a boundary layer. This layer is a thin region at the right side of the domain where the solution has a steep gradient. As the solution gets steeper, the classic methods on the uniform mesh are unable to solve SPDDEs without using an inadequately small mesh size, which is not practicable. This motivated various researchers to develop ε-uniform non-classical numerical methods. When the delay parameter τ is smaller than ε, the given delay differential equation is reduced by means of Taylor series expansion, but if the delay parameter is higher, this approach does not work [10]. In such case, during developing the scheme the point t − τ must coincide with a mesh point. In this study, we have a numerical scheme for the problem with the large delay.

In the last few years, several numerical approaches have been developed for the solution of time dependent singularly perturbed differential equations [8]. However, in the case of these equations with time delay, few numerical methods have been developed, and studies are still at an early stage [8]. Some of them, which treat time dependent singularly perturbed delay parabolic convection-diffusion initial boundary value problems are listed here. Das and Natesan [11] had proposed to a numeric scheme consisting of an implicit-Euler scheme in the temporal direction, accompanied by a hybrid scheme in the space direction. Negero and Duressa [12] developed a second-order ε-uniform convergent scheme on a uniform mesh. They discretize time and spatial derivatives by the implicte Euler rule and Micken's non-standard method with extrapolation, respectively. In the same year, Negero and Duressa [13] designed an efficient numerical approach which is uniformly convergent of second order of convergent. They discretize time and spatial derivatives, respectively, by Crank-Nicolson and the exponentially fitted spline method. Woldaregay et al. [14] developed a novel numerical scheme using by an exponentially fitted operator which is ε-uniform convergent of linear order of convergent. Kumar et al. [15] developed a graded mesh refinement approach. The approach shows parameter uniform convergent linear order. Abdelhakem and Youssri [16] proposed two spectral Legendre's derivative algorithms for Lane-Emden, Bratu equations, and singular perturbed problems. They have shown the numerical schemes are stable, convergent, and accurate. Abd-Elhameed et al. [17] suggested and analyzed a new operational matrix method based on shifted Legendre polynomials for obtaining numerical spectral solutions of linear and non-linear second-order boundary value problems. The authors showed that the method has the following advantages. The method can be applied for both linear and non-linear second-order boundary value problems including some important singular perturbed equations and also a Bratu-type equation. This method computes highly accurate approximate solutions. Authors [18–24] also developed ε-uniform convergent numerical methods for singularly perturbed delay partial differential equation.

Nowadays, the use of the spline based approach has become very popular among different numerical methods to solve SPDDEs. Daba and Duressa [25] proposed a uniform convergent numerical scheme for singularly perturbed parabolic convection- diffusion equation with a small delay and advance parameter in the spatial variable of reaction term using an extended cubic B-spline method. They also suggested uniformly convergent numerical scheme for this problem using cubic B-spline [26] on a uniform mesh. Kumar and Kadalbajoo [27] proposed the parameter-uniform numerical method for the problem using cubic B-spline on a Shishkin mesh. Kumar [8] and Negero and Duressa [28] developed a parameter uniform convergent method to solve time dependent singularly perturbed delay parabolic convection-diffusion initial boundary value problems using the cubic B-spline collocation method on a piecewise uniform Shishkine mesh and uniform mesh, respectively. The extended cubic B-spline collocation method is a generalization of the cubic B-spline collocation method. It introduces a free parameter to allow the cubic B-spline's shape to alter while maintaining the continuity in the order of three. Kumar and Kumari [29] designed ε-uniform convergent numerical methods using an extended cubic B-spline collocation method for singularly perturbed delay parabolic convection-diffusion initial boundary value problems on a piecewise uniform Shishkine mesh.

The main contribution of this study is to construct a parameter uniform numerical method for singularly perturbed delay parabolic convection-diffusion problems. The method consists of the implicit Euler rule for temporal discretization and the extended cubic B-spline collocation method for spatial discretization on a uniform mesh using artificial viscosity. In this method, we use artificial viscosity σ(x, ε) to substitute the perturbation parameter ε that affects the highest derivative.

The rest of this article is organized as follows: in Section 2, we discuss the continuous problem and show the boundedness of the exact solution. The numerical scheme of the problem is discussed in Section 3. In this section, the implicit Euler method, the extended cubic B-spline collocation method, and the derivation of artificial viscosity are briefly discussed. Section 4 describes the convergence analysis of the proposed method. In Section 5, numerical illustrations are given to confirm the theoretical investigation. Lastly, Section 6 gives conclusion of the article.

Notation

C, C₁, C₂: A generic positive constants independent of ε and the mesh parameters N and M

C^k(D): The set of k times continuously differentiable function on domain D

∥f∥_D ${sup}_{\left(x,t\right)\in D}\left|f\left(x,t\right)\right|$

2. Problem statement

Let Ω_x = (0, 1), Ω_t* = [−τ, 0] and Ω_t = (0, T] for some fixed positive time T. We define D = Ω_x × Ω_t and ∂D = D₀ ∪ D_b ∪ D₁, where $D_0=\left\{\left(0,t\right),t\in \overline{{\Omega }_t}\right\}$, $D_1=\left\{\left(1,t\right),t\in \overline{{\Omega }_t}\right\}$, and $D_b=\overline{{\Omega }_x}\times {\Omega }_{t*}.$ We consider the following singularly perturbed delay parabolic initial boundary value problem (IBVP):

$\begin{array}{ll}\left(\frac{\partial }{\partial t}+{𝔏}_{\epsilon }\right)u\left(x,t\right)=-c\left(x,t\right)u\left(x,t-\tau \right)+f\left(x,t\right)\left(x,t\right)\in D& (1)\end{array}$

with

$\begin{array}{ll}u\left(0,t\right)={\psi }_0\left(t\right),u\left(1,t\right)={\psi }_1\left(t\right)t\in \overline{{\Omega }_t}& (2)\end{array}$

and

$\begin{array}{ll}u\left(x,t\right)={\psi }_b\left(x,t\right),\left(x,t\right)\in D_b& (3)\end{array}$

where τ and 0 < ε ≪ 1 are delay and singular perturbation parameter, respectively. The differential operator 𝔏_ε in Eq. (1) is defined as

${𝔏}_{\epsilon }u\left(x,t\right)=-\epsilon \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+a\left(x\right)\frac{\partial u\left(x,t\right)}{\partial x}+b\left(x,t\right)u\left(x,t\right).$

The functions a(x), b(x, t), c(x, t), f(x, t), ψ₀(t), ψ₁(t), and ψ_b(x, t) are assumed to be sufficiently smooth, bounded, and independent of ε. It is also assumed that

$a\left(x\right)\ge \alpha >0,b\left(x,t\right)\ge \beta >0,c\left(x,t\right)\ge \gamma >0,\left(x,t\right)\in \overline{D}$

and terminal T, satisfy T = kτ for some positive integer k. We also assume that the data satisfy the following compatibility conditions:

$\begin{array}{l}\begin{array}{c}{\psi }_b\left(0,0\right)={\psi }_0\left(0\right),\hfill \\ {\psi }_b\left(1,0\right)={\psi }_1\left(0\right),\hfill \\ \frac{d{\psi }_0}{dt}{|}_{t=0}-\epsilon \frac{{\partial }^2{\psi }_b}{\partial x^2}{|}_{\left(0,0\right)}+a\left(0\right)\frac{\partial {\psi }_b}{\partial x}{|}_{\left(0,0\right)}\hfill \\ +b\left(0,0\right){\psi }_b\left(0,0\right)=-c\left(0,0\right){\psi }_b\left(0,-\tau \right)+f\left(0,0\right),\hfill \\ \frac{d{\psi }_1}{dt}{|}_{t=0}-\epsilon \frac{{\partial }^2{\psi }_b}{\partial x^2}{|}_{\left(1,0\right)}+a\left(1\right)\frac{\partial {\psi }_b}{\partial x}{|}_{\left(1,0\right)}\hfill \\ +b\left(1,0\right){\psi }_b\left(1,0\right)=-c\left(1,0\right){\psi }_b\left(1,-\tau \right)+f\left(1,0\right),\hfill \end{array}\hfill \end{array}$

Under the above assumptions and conditions, the solution of IBVP Eq. (1) is unique which has a parabolic boundary layer of width $\mathcal{O}\left(\epsilon \right)$ along x = 1. Compatibility conditions [30] are relationships between the data of the problem and the differential operator that ensure that derivatives of u(x, t) up to a desired order are continuous on the closed domain D. They do not result from the problem's singularly perturbed nature because they only appear at corners.

Lemma 2.1 (Continuous Maximum Principle). Suppose the function $\phi \left(x,t\right)\in C^{2,1}\left(\overline{D}\right)$ satisfies $\left(\frac{\partial }{\partial t}+{𝔏}_{\epsilon }\right)\phi \left(x,t\right)\ge 0,\forall \left(x,t\right)\in D$ and φ(x, t) ≥ 0, ∀(x, t) ∈ ∂D, then $\phi \left(x,t\right)\ge 0,\forall \left(x,t\right)\in \overline{D}.$

Proof. Let $\left(\tilde{x},\tilde{t}\right)\in \overline{D}$, such that $\phi \left(\tilde{x},\tilde{t}\right)={min}_{\left(x,t\right)\in \overline{D}}\phi \left(x,t\right)$ and suppose that $\phi \left(\tilde{x},\tilde{t}\right)<0.$ Clearly $\left(\tilde{x},\tilde{t}\right)\notin \partial D$. Then, it follows from calculus that $\frac{\partial \phi \left(\tilde{x},\tilde{t}\right)}{\partial x}=0,\frac{\partial \phi \left(\tilde{x},\tilde{t}\right)}{\partial t}=0$ and $\frac{{\partial }^2\phi \left(\tilde{x},\tilde{t}\right)}{\partial x^2}\ge 0$. Therefore, from Eq (1), we have

$\begin{array}{l}\left(\frac{\partial }{\partial t}+{𝔏}_{\epsilon }\right)\phi \left(\tilde{x},\tilde{t}\right)=\frac{\partial \phi \left(\tilde{x},\tilde{t}\right)}{\partial t}-\epsilon \frac{{\partial }^2\phi \left(\tilde{x},\tilde{t}\right)}{\partial x^2}+a\left(\tilde{x}\right)\frac{\partial \phi \left(\tilde{x},\tilde{t}\right)}{\partial x}\\ +b\left(\tilde{x},\tilde{t}\right)\phi \left(\tilde{x},\tilde{t}\right)\le 0,\end{array}$

which is a contradiction to the assumption made. Thus, $\phi \left(\tilde{x},\tilde{t}\right)\ge 0$ which leads to $\phi \left(x,t\right)\ge 0,\forall \left(x,t\right)\in \overline{D}$.     □

Lemma 2.2. [11] The solution u(x, t) of the IBVP Eq. (1) satisfies the estimate

$\begin{array}{ll}|u\left(x,t\right)-{\psi }_b\left(x,0\right)|\le Ct,\forall \left(x,t\right)\in \overline{D}& (4)\end{array}$

Lemma 2.3. The solution u(x, t) of the IBVP Eq. (1) satisfies the estimate

$\begin{array}{ll}\left|u\left(x,t\right)\right|\le C,\forall \left(x,t\right)\in \overline{D}.& (5)\end{array}$

Proof. Since t ∈ (0, T], by Lemma 2.2, we get |u(x, t) − ψ_b(x, 0)| ≤ Ct ≤ CT. It follows that

$\begin{array}{l}\left|u\left(x,t\right)\right|=|u\left(x,t\right)-{\psi }_b\left(x,0\right)+{\psi }_b\left(x,0\right)|\\ \le |u\left(x,t\right)-{\psi }_b\left(x,0\right)|+\left|{\psi }_b\left(x,0\right)\right|,\\ \le CT+\left|{\psi }_b\left(x,0\right)\right|.\end{array}$

Since $\left|{\psi }_b\left(x,0\right)\right|\in C^{2,1}\left(\overline{D}\right)$, CT+|ψ_b(x, 0)| is bounded by some positive constant C, $\left|u\left(x,t\right)\right|\le C,\forall \left(x,t\right)\in \overline{D}$.     □

Lemma 2.4 (Uniform stability estimate). Let u(x, t) be solution of the IBVP Eq. (1). Then, we obtain the bound:

$\begin{array}{l}\left|u\left(x,t\right)\right|\le {\beta }^{-1}\Vert \left(\frac{\partial }{\partial t}+{𝔏}_{\epsilon }\right)u\left(x,t\right)\Vert \\ +max\left\{\left|{\psi }_0\left(t\right)\right|,\left|{\psi }_1\left(t\right)\right|,\left|{\psi }_b\left(x,t\right)\right|\right\}.& (6)\end{array}$

Proof. Let define barrier functions

${\text{Ψ}}^{\pm }\left(x,t\right)={\beta }^{-1}\left|\right|\left(\frac{\partial }{\partial t}+{𝔏}_{\epsilon }\right)u\left(x,t\right)\left|\right|+max\left\{\left|{\psi }_0\left(t\right)\right|,\left|{\psi }_1\left(t\right)\right|,\left|{\psi }_b\left(x,t\right)\right|\right\}\pm u\left(x,t\right)$. By applying the continuous maximum principle Lemma 2.1, we obtain the required result.     □

Theorem 2.1. The solution u(x, t) of Eq (1) and its derivatives satisfy the following bounds :

$\begin{array}{ll}\left|\frac{{\partial }^{i+j}u\left(x,t\right)}{\partial x^j\partial t^i}\right|\le C\left(1+{\epsilon }^{-j}exp\left(-\alpha \left(1-x\right)/\epsilon \right)\right),\text{for all},\left(x,t\right)\in \overline{D},& (7)\end{array}$

where i and j are non-negative integers such that 0 ≤ i + j ≤ 5.

Proof. See details of the proof in Mbroh et al. [31].     □

3. Numerical scheme

This section describes the semi-discretization and extended cubic B-spline method on the uniform mesh by introducing artificial viscosity.

3.1. The time semi-discretization

Here, we use Rothe's technique [32, 33] to discretize the time variable of equation Eq. (1) by means of the Euler implicit rule. We make the argument t − τ as a nodal point for the uniform mesh of step size Δt. Let ${\Omega }_t^M=\left\{t_i=i\Delta t=i\left(T/M\right),i=0,1,\dots ,M,t_M=T\right\}$ and ${\Omega }_t^m=\left\{-t_i=i\Delta t,i=0,1,\dots ,m,t_m=\tau \right\}$ for the intervals [0, T] and [−τ, 0], respectively. M and m are the number of mesh elements to their respective intervals. After time discretization, we obtain the following semi-discretized problem

$\begin{array}{ll}\{\begin{array}{cc}{𝔏}_{\epsilon }^M{U}^{i+1}\left(x\right)={\Phi }^i\left(x\right),\hfill & x\in {\Omega }_x,0\le i\le M-1\\ U^{i+1}\left(0\right)={\psi }_0\left(t_{i+1}\right),\hfill & 0\le i\le M-1\\ U^{i+1}\left(1\right)={\psi }_1\left(t_{i+1}\right),\hfill & 0\le i\le M-1\\ U^{-i}\left(x,t\right)={\psi }_b\left(x,-t_i\right),\hfill & x\in \overline{{\Omega }_x},0\le i\le m\end{array}& (8)\end{array}$

where ${𝔏}_{\epsilon }^M{U}^{i+1}\left(x\right)\equiv -\epsilon \frac{d^2{U}^{i+1}\left(x\right)}{d{x}^2}+a\left(x\right)\frac{d{U}^{i+1}\left(x\right)}{dx}+d^{i+1}{U}^{i+1}\left(x\right)$, $d^{i+1}\left(x\right)=\frac{1}{\Delta t}+b^{i+1}\left(x\right)$, ${\Phi }^i\left(x\right)=\frac{1}{\Delta t}{U}^i\left(x\right)-c^{i+1}\left(x\right)U^{i+1-m}\left(x\right)+f^{i+1}\left(x\right)$.

Uⁱ⁺¹(x) is the approximation of the exact solution u(x, t_i+1) at (i + 1)^th time level.

Lemma 3.1 (Semi-discrete maximum principle). Let yⁱ⁺¹(x) be a sufficiently smooth function on $\overline{{\Omega }_x}$ such that yⁱ⁺¹(0) ≥ 0, yⁱ⁺¹(1) ≥ 0. Then, ${𝔏}_{\epsilon }^M{y}^{i+1}\left(x\right)\ge 0,\forall x\in {\Omega }_x$ implies yⁱ⁺¹(x) ≥ 0, and $\forall x\in \overline{{\Omega }_x}$.

Proof. Let $x^{*}\in \overline{{\Omega }_x}$ such that $\underset{x\in \overline{{\Omega }_x}}{min}{y}^{i+1}\left(x\right)=y^{i+1}\left(x^{*}\right)$ and assume yⁱ⁺¹(x^*) < 0. It is clear that x^* ∉ {0, 1}. From property of calculus, we have $\frac{d{y}^{i+1}\left(x^{*}\right)}{dx}=0$ and $\frac{d^2{y}^{i+1}\left(x^{*}\right)}{d{x}^2}\ge 0$. This yield ${𝔏}_{\epsilon }^M{y}^{i+1}\left(x^{*}\right)<0$ which contradicts to the hypothesis made ${𝔏}_{\epsilon }^M{y}^{i+1}\left(x\right)\ge 0$. Therefore, we conclude that $y^{i+1}\left(x\right)\ge 0,\forall x\in \overline{{\Omega }_x}$.     □

In the temporal semi-discretization, $e_{i+1}=u\left(x,t_{i+1}\right)-{\hat{U}}^{i+1}\left(x\right)$ defines the local truncation error, where ${\hat{U}}^{i+1}\left(x\right)$ is the solution obtained after one step of the semi-discrete scheme taking the exact solution U(x, t_i) instead of Uⁱ(x) as the starting data. For each time step, local error estimate contribute to the global error in temporal discretization which is defined at the instant t_i as $E_i=u\left(x,t_i\right)-U^i\left(x\right)$.

Lemma 3.2 (Local error estimate). The local error corresponding to the semi-discretized problem Eq. (8) satisfies

$\Vert e_{i+1}\Vert \le C{\left(\Delta t\right)}^2$

Proof. Taylor's series expansion on u(x, t) gives

$\begin{array}{ll}u\left(x,t_{i+1}\right)-\Delta t\frac{\partial u\left(x,t_{i+1}\right)}{\partial t}=u\left(x,t_i\right)+\mathcal{O}{\left(\Delta t\right)}^2& (9)\end{array}$

Substituting Eq. (1) in to Eq. (9), we get

$\begin{array}{l}u(x,t_{i+1})-\Delta t[\epsilon \frac{{\partial }^{2}u(x,t_{i+1})}{\partial x^2}-a(x)\frac{\partial u(x,t_{i+1})}{\partial t}\\ -b(x,t_{i+1})u(x,t_{i+1})-c(x,t_{i+1})u(x,t_{i+1-m})+f(x,t_{i+1})]\\ =u(x,t_i)+\mathcal{O}{(\Delta t)}^2& (10)\end{array}$

In terms of the differential operator, we get

$\begin{array}{l}(1+\Delta t{L}_{\epsilon })u(x,t_{i+1})=u(x,t_i)+\Delta t(-c(x,t_{i+1})u(x,t_{i+1-m})\\ +f(x,t_{i+1}))+\mathcal{O}{(\Delta t)}^2& (11)\end{array}$

Thus, ${\hat{U}}^{i+1}\left(x\right)$ satisfies

$\begin{array}{ll}\left(1+\Delta t{𝔏}_{\epsilon }\right){\hat{U}}^{i+1}\left(x\right)={\hat{U}}^i\left(x\right)+\Delta t\left(-c\left(x,t_{i+1}\right){\hat{U}}^{i+1-m}\left(x\right)+f\left(x,t_{i+1}\right)\right).\hfill & (12)\end{array}$

From Eq. (11) and Eq. (12), the local error satisfy the following boundary value problem:

$\begin{array}{ll}\{\begin{array}{c}\left(1+\Delta t{𝔏}_{\epsilon }\right)e_{i+1}=\mathcal{O}\left({\left(\Delta t\right)}^2\right)\hfill \\ e_{i+1}\left(0\right)=0=e_{i+1}\left(1\right).\hfill \end{array}\hfill & (13)\end{array}$

An application of Lemma 3.1 on the operator (1 + Δt𝔏_ε) gives $\left|\right|e_{i+1}\left|\right|\le C{\left(\Delta t\right)}^2$, which completes the proof.     □

Lemma 3.3 (Global error estimate). The global error estimate at t_i+1 satisfies

‖E_i+1‖ ≤ CΔt, iΔt ≤ T.

Proof. By definition, we have

$\begin{array}{l}\Vert E_{i+1}\Vert =\Vert {\displaystyle \sum _{s=1}^i}{e}_s\Vert \\ \le \Vert e_1\Vert +\Vert e_2\Vert +\cdots +\Vert e_i\Vert \\ \le iC\left({\left(\Delta t\right)}^2\right),\text{Lemma }3.2\\ =C\left(i\Delta t\right)\Delta t,i\Delta t\le T\\ \le C\Delta t,\end{array}$

where C is a positive constant independent of ε and Δt.     □

Lemma 3.4. The solution of Eq. (8) satisfies

$\begin{array}{l}\left|\frac{d^n{U}^{i+1}\left(x\right)}{d{x}^n}\right|\le C\left(1+{\epsilon }^{-n}exp\left(-\frac{\alpha \left(1-x\right)}{\epsilon }\right)\right),\\ x\in \overline{{\Omega }_x},n=0,1,2,3,4.\end{array}$

Proof. See the proof in Kellogg and Tsan [34] and Clavero et al. [35].     □

3.2. The spatial discretization

Here, we apply the extended cubic B-spline collocation for the problem Eq. (8). Artificial viscosity shall be introduced to take into account the exponential properties of exact solution on the uniform mesh. Therefore, the perturbation parameter ε, which affects the highest derivative, is replaced by an artificial viscosity σ(x, ε). We rewrite the problem Eq. (8) as

$\begin{array}{ll}\{\begin{array}{cc}{𝔏}_{\sigma }^M{U}^{i+1}\left(x\right)={\Phi }^i\left(x\right),\hfill & x\in {\Omega }_x,0\le i\le M-1\\ U^{i+1}\left(0\right)={\psi }_0\left(t_{i+1}\right),\hfill & 0\le i\le M-1\\ U^{i+1}\left(1\right)={\psi }_1\left(t_{i+1}\right),\hfill & 0\le i\le M-1\\ U^{-i}\left(x,t\right)={\psi }_b\left(x,-t_i\right),\hfill & x\in \overline{{\Omega }_x},0\le i\le m,\end{array}\hfill & (14)\end{array}$

where ${𝔏}_{\sigma }^M{U}^{i+1}\left(x\right)\equiv -\sigma \left(x,\epsilon \right)\frac{d^2{U}^{i+1}\left(x\right)}{d{x}^2}+a\left(x\right)\frac{d{U}^{i+1}\left(x\right)}{dx}+d^{i+1}{U}^{i+1}\left(x\right)$, $d^{i+1}\left(x\right)=\frac{1}{\Delta t}+b^{i+1}\left(x\right)$, ${\Phi }^i\left(x\right)=\frac{1}{\Delta t}{U}^i\left(x\right)-c^{i+1}\left(x\right)U^{i+1-m}\left(x\right)+f^{i+1}\left(x\right)$. Now, Uⁱ⁺¹(x) is the approximation of the exact solution u(x, t_i+1) at the (i + 1)^th time level after introducing artificial viscosity. The properties of all data in Eq. (8) are retained in Eq. (14).

3.2.1. Extended cubic B-spline method

The interval [0, 1] is divided such that knots are equally distributed as ${\Omega }_x^N={\left\{x_i\right\}}_0^N$ and mesh spacing $h=\frac{1}{N}$. Let η ∈ ℝ, then the blending function of degree 4 of extended cubic B-spline Q_j has the following form [36]:

$\begin{array}{ll}{Q}_j\left(x,\eta \right)=\frac{1}{24h^4}\{\begin{array}{cc}4h\left(1-\eta \right){\left(x-x_{j-2}\right)}^3+3\eta {\left(x-x_{j-2}\right)}^4,\hfill & x\in \left[x_{j-2},x_{j-1}\right)\\ h^4\left(4-\eta \right)+12h^3\left(x-x_{j-1}\right)+6h^2\left(2+\eta \right){\left(x-x_{j-1}\right)}^2\hfill \\ -12h{\left(x-x_{j-1}\right)}^3-3\eta {\left(x-x_{j-1}\right)}^4,\hfill & x\in \left[x_{j-1},x_j\right)\\ h^4\left(4-\eta \right)+12h^3\left(x_{j+1}-x\right)+6h^2\left(2+\eta \right){\left(x_{j+1}-x\right)}^2\hfill \\ -12h{\left(x_{j+1}-x\right)}^3-3\eta {\left(x_{j+1}-x\right)}^4,\hfill & x\in \left[x_j,x_{j+1}\right)\\ 4h\left(1-\eta \right){\left(x_{j+2}-x\right)}^3+3\eta {\left(x_{j+2}-x\right)}^4,\hfill & x\in \left[x_{j+1},x_{j+2}\right)\\ 0,\hfill & \text{Otherwise}\end{array}\hfill & (15)\end{array}$

where η is a free parameter which is used to obtain different form of extended cubic B-spline functions and satisfy [37] the property, −8 ≤ η ≤ 1. When η = 0, Q_j(x, η) degenerates into exactly cubic B-spline functions. The extended cubic B-spline is the generalization of the B-spline. The sequence $\Lambda ={\left\{Q_j\left(x,\eta \right)\right\}}_{-1}^{N+1}$ forms a basis for the functions defined over the interval [0, 1] and $Q_j\left(x,\eta \right)\in C^2\left({\Omega }_x\right)$. The value of $Q_j\left(x,\eta \right),Q_j^{{}^{\prime }}\left(x,\eta \right)$, and $Q_j^{{}^{″}}\left(x,\eta \right)$ at the knots x_js computed from Eq. (15) are shown Table 1.

TABLE 1

Table 1. Values of $Q_j,Q_j^{{}^{\prime }},Q_j^{{}^{″}}$ at knots.

Now, suppose that the approximate solution $S\left(x,\eta \right)\in {\hat{Q}}_3\left({\Omega }_x^N\right)=span\Lambda$ to the exact solution u(x, t_i+1) at (i + 1)^th time level given by

$\begin{array}{ll}S\left(x,\eta \right)={\displaystyle \sum _{j=-1}^{N+1}}{\delta }_j{Q}_j\left(x,\eta \right),& (16)\end{array}$

where δ_j are the parameters to be determined from the application of the collocation method, initial and boundary conditions. Using Table 1, evaluation of Eq. (16) and its first and second derivatives at knots x_j yield

$\begin{array}{ll}\begin{array}{c}S\left(x_j,\eta \right)=\frac{4-\eta }{24}{\delta }_{j-1}+\frac{8+\eta }{12}{\delta }_j+\frac{4-\eta }{24}{\delta }_{j+1}\hfill \\ S^{{}^{\prime }}\left(x_j,\eta \right)=-\frac{1}{2h}\left({\delta }_{j-1}-{\delta }_{j+1}\right)\hfill \\ S^{{}^{″}}\left(x_j,\eta \right)=\frac{2+\eta }{2h^2}\left({\delta }_{j-1}-2{\delta }_j+{\delta }_{j+1}\right).\hfill \end{array}\hfill & (17)\end{array}$

Using Eq. (17) at x_j in Eq. (14), we get the system of N + 1 linear equations in N + 3 unknown as

$\begin{array}{c}-{\sigma }_j\left[\frac{2+\eta }{2h^2}\left({\delta }_{j-1}-2{\delta }_j+{\delta }_{j+1}\right)\right]+a_j\left[-\frac{1}{2h}\left({\delta }_{j-1}-{\delta }_{j+1}\right)\right]\hfill \\ +d_j^{i+1}\left[\frac{4-\eta }{24}{\delta }_{j-1}+\frac{8+\eta }{12}{\delta }_j+\frac{4-\eta }{24}{\delta }_{j+1}\right]={\Phi }^i\left(x_j\right).\hfill & (18)\end{array}$

This gives

$\begin{array}{ll}{R}_j^{-}{\delta }_{j-1}+R_j^0{\delta }_j+R_j^{+}{\delta }_{j+1}=h^2{\Phi }_j^i,0\le j\le N,\hfill & (19)\end{array}$

where

$\begin{array}{l}\{\begin{array}{c}{R}_j^{-}=-{\sigma }_j\left(\frac{2+\eta }{2}\right)-\frac{1}{2}h{a}_j+\left(\frac{4-\eta }{24}\right)h^2{d}_j^{i+1}\hfill \\ R_j^0={\sigma }_j\left(2+\eta \right)+\left(\frac{8+\eta }{12}\right)h^2{d}_j^{i+1}\hfill \\ R_j^{+}=-{\sigma }_j\left(\frac{2+\eta }{2}\right)+\frac{1}{2}h{a}_j+\left(\frac{4-\eta }{24}\right)h^2{d}_j^{i+1}\hfill \end{array}\hfill \end{array}$

σ_j = σ(x_j, ε),      a_j = a(x_j),

3.2.2. Design of the artificial viscosity

The two-variable method [38] for boundary value problem Eq. (14) with a(x) > 0 expresses the solution Uⁱ⁺¹(x) as

$\begin{array}{l}{U}^{i+1}\left(x\right)=U_0^{i+1}\left(x\right)+\frac{a\left(1\right)}{a\left(x\right)}\left({\psi }_1\left(t_{i+1}\right)-U_0^{i+1}\left(1\right)\right)\hfill \\ exp\left(-{\displaystyle {\int }_x^1}\frac{a^2\left(s\right)+\epsilon d^{i+1}\left(s\right)}{\epsilon a\left(s\right)}ds\right)+\mathcal{O}\left(\epsilon \right),& (20)\end{array}$

where $U_0^{i+1}\left(x\right)$ is solution of the reduced problem

$a\left(x\right)\frac{d{U}^{i+1}\left(x\right)}{dx}+d^{i+1}\left(x\right)U^{i+1}\left(x\right)={\Phi }^i\left(x\right),U_0^{i+1}\left(0\right)=U_0\left(0\right).$

With the Taylor series expanding a(x) and dⁱ⁺¹(x) about “1," and keeping the first term, we obtain

$\begin{array}{l}{U}^{i+1}\left(x\right)=U_0^{i+1}\left(x\right)+\left({\psi }_1\left(t_{i+1}\right)-U_0^{i+1}\left(1\right)\right)\\ exp\left(-\frac{\left(a^2\left(1\right)+\epsilon d^{i+1}\left(1\right)\right)\left(1-x\right)}{\epsilon a\left(1\right)}\right)+\mathcal{O}\left(\epsilon \right).& (21)\end{array}$

Eqs (19) and (21) at the nodal point x_j are given by

$\begin{array}{l}\left[-\left(\frac{2+\eta }{2}\right)\frac{{\sigma }_j}{h}-\frac{a(x_j)}{2}+\left(\frac{4-\eta }{24}\right)h{d}_j^{i+1}\right]{\delta }_{j-1}+[(2+\eta )\frac{{\sigma }_j}{h}\\ +\left(\frac{8+\eta }{12}\right)h{d}_j^{i+1}]{\delta }_j\\ +\left[-\left(\frac{2+\eta }{2}\right)\frac{{\sigma }_j}{h}+\frac{a(x_j)}{2}+\left(\frac{4-\eta }{24}\right)h{d}_j^{i+1}\right]{\delta }_{j+1}=h{\Phi }_j^i\end{array}$

and

$\begin{array}{c}\frac{4-\eta }{24}{\delta }_{j-1}+\frac{8+\eta }{12}{\delta }_j+\frac{4-\eta }{24}{\delta }_{j+1}=U_0^{i+1}\left(x\right)\hfill \\ +\left({\psi }_1\left(t_{i+1}\right)-U_0^{i+1}\left(1\right)\right)exp\left(-\frac{\left(a^2\left(1\right)+\epsilon d^{i+1}\left(1\right)\right)\left(1-x\right)}{\epsilon a\left(1\right)}ds\right)+\mathcal{O}\left(\epsilon \right).\hfill \end{array}$

In th limiting case as h → 0, we get

$R=\left[\begin{array}{ccccccc}{R}_0^0-2\left(\frac{8+\eta }{4-\eta }\right)R_0^{-}& R_0^{+}-R_0^{-}& & & & & \\ R_1^{-}& R_1^0& R_1^{+}& & & & \\ & \ddots & \ddots & \ddots & & & \\ & & R_j^{-}& R_j^0& R_j^{+}& & \\ & & & \ddots & \ddots & \ddots & \\ & & & & R_{N-1}^{-}& R_{N-1}^0& R_{N-1}^{+}\\ & & & & & R_N^{-}-R_N^{+}& R_N^0-2\left(\frac{8+\eta }{4-\eta }\right)R_N^{+}\end{array}\right]$

$\begin{array}{ll}{\displaystyle \underset{h\to 0}{lim}}\frac{{\sigma }_j}{h}=\frac{a\left(0\right)}{2+\eta }\left(\frac{{\delta }_{j-1}-{\delta }_{j+1}}{{\delta }_{j-1}-2{\delta }_j+{\delta }_{j+1}}\right)& (22)\end{array}$

and

$\begin{array}{c}\frac{4-\eta }{24}{\delta }_{j-1}+\frac{8+\eta }{12}{\delta }_j+\frac{4-\eta }{24}{\delta }_{j+1}=U_0^{i+1}\left(jh\right)\hfill \\ +\left({\psi }_1\left(t_{i+1}\right)-U_0^{i+1}\left(1\right)\right)exp\left(-\frac{\left(a^2\left(1\right)+\epsilon d^{i+1}\left(1\right)\right)\left(1/\epsilon -i\rho \right)}{a\left(1\right)}\right),\hfill & (23)\end{array}$

where ρ = h/ε. Similarly, computing the values of $\underset{h\to 0}{\lim }\frac{{\sigma }_j}{h}$ at x_j−1, x_j, x_j+1 and adding in the proportion $\frac{4-\eta }{24},\frac{8+\eta }{12},\frac{4-\eta }{24}$, respectively, and removing δ_j's using Eq. (23), we define

$\begin{array}{ll}{\sigma }_j=\epsilon \frac{a\left(x_j\right)\rho }{2+\eta }coth\left(\frac{a\left(x_j\right)\rho }{2}\right)& (24)\end{array}$

Since, $cothx-\frac{1}{x}=\frac{x}{3}+\mathcal{O}\left(x^3\right)\mathrm{\text{as}}x\to 0$ and $cothx-\frac{1}{x}=1+\mathcal{O}\left(\frac{1}{x}\right)\mathrm{\text{as}}x\to \infty$. We have

$C_1\frac{x^2}{x+1}\le xcothx-1\le C_2\frac{x^2}{x+1},\epsilon \frac{{\left(h/\epsilon \right)}^2}{h/\epsilon +1}=\frac{h^2}{h+\epsilon }$

Thus, $|xcothx-1|\le C\frac{x^2}{x+1}$, which implies

$\begin{array}{ll}|{\sigma }_j-\epsilon |\le C\frac{h}{h+\epsilon }& (25)\end{array}$

For 0 ≤ j ≤ N, Eq. (19) is a system of (N + 1) linear equation in (N + 3) unknown δ₋₁, δ₀, δ₁, …, δ_N+1. Now, by applying boundary conditions, at i = 0 and i = N, for Eq. (19) and first equation of Eq. (17), we can eliminate δ₋₁ and δ_N+1. Therefore, we get

$\begin{array}{l}\left(R_0^0-2\left(\frac{8+\eta }{4-\eta }\right)R_0^{-}\right){\delta }_0+\left(R_0^{+}-R_0^{-}\right){\delta }_1\\ =h^2{\Phi }_0^i-\left(\frac{24}{4-\eta }\right)R_0^{-}{\psi }_0\left(t_{i+1}\right),& (26)\end{array}$

$\begin{array}{l}\left(R_N^{-}-R_N^{+}\right){\delta }_{N-1}+\left(R_N^0-2\left(\frac{8+\eta }{4-\eta }\right)R_N^{+}\right){\delta }_N\\ =h^2{\Phi }_N^i-\left(\frac{24}{4-\eta }\right)R_N^{+}{\psi }_1\left(t_{i+1}\right).& (27)\end{array}$

Finally, in (N + 1) unknowns δ₀, δ₁, …, δ_N, we find a system of (N + 1) equations with a matrix form as

$\begin{array}{ll}R\delta =\Phi & (28)\end{array}$

where

$\delta =\left[\begin{array}{c}{\delta }_0\\ {\delta }_1\\ ⋮\\ ⋮\\ ⋮\\ {\delta }_{N-1}\\ {\delta }_N\end{array}\right],\Phi =\left[\begin{array}{c}{h}^2{\Phi }_0^i-\left(\frac{24}{4-\eta }\right)R_0^{-}{\psi }_0\left(t_{i+1}\right)\\ h^2{\Phi }_1^i\\ ⋮\\ ⋮\\ ⋮\\ h^2{\Phi }_{N-1}^i\\ h^2{\Phi }_N^i-\left(\frac{24}{4-\eta }\right)R_N^{+}{\psi }_1\left(t_{i+1}\right)\end{array}\right]$

As h → 0, R is diagonally dominant tridiagonal matrix, non-singular. Therefore, we can compute the values of δ_j's, then substituted in Eq. (16), to obtain the approximate solution to Eq. (14).

4. Parameter-uniform convergence analysis

In this section, we establish the parameter-uniform convergence of the extended cubic B-spline collocation method. The following lemma will be used in the convergence analysis.

Lemma 4.1. The extended cubic B-spline set Λ = {Q₋₁, Q₀, Q₁, …, Q_N+1} defined in Eq. (15), satisfy the inequality

$\begin{array}{l}{\displaystyle \sum _{j=-1}^{N+1}}\left|Q_j\left(x,\eta \right)\right|\le 1.75,0\le x\le 1.\end{array}$

Proof. We know that

$\left|{\displaystyle \sum _{j=-1}^{N+1}}{Q}_j\left(x_i,\eta \right)\right|\le {\displaystyle \sum _{j=-1}^{N+1}}|Q_j\left(x_i,\eta \right)|.$

Let x = x_i be a nodal point. Then, by the definition of Q_j(x, η), we have

$\begin{array}{l}{\displaystyle \sum _{j=-1}^{N+1}}\left|Q_j\left(x_i,\eta \right)\right|=\left|Q_{i-1}\left(x_i,\eta \right)\right|+\left|Q_i\left(x_i,\eta \right)\right|+\left|Q_{i+1}\left(x_i,\eta \right)\right|\\ =\frac{4-\eta }{24}+\frac{8+\eta }{12}+\frac{4-\eta }{24}=1<1.75.\end{array}$

Now, for x_j−1 < x < x_j, from Table 1, we get

$\begin{array}{l}{Q}_j\left(x,\eta \right)<\frac{8+\eta }{12},Q_{j-1}\left(x,\eta \right)=\frac{8+\eta }{12}\\ Q_{j+1}\left(x,\eta \right)<\frac{4-\eta }{24},Q_{j-2}\left(x,\eta \right)=\frac{4-\eta }{24}\end{array}$

Thus,

$\begin{array}{l}{\displaystyle \sum _{j=-1}^{N+1}}\left|Q_j\left(x,\eta \right)\right|=\left|Q_{j-1}\left(x,\eta \right)\right|+\left|Q_j\left(x,\eta \right)\right|+\left|Q_{j+1}\left(x,\eta \right)\right|\\ +\left|Q_{j-2}\left(x,\eta \right)\right|\le \frac{20+\eta }{12}\end{array}$

Since, $-8\le \eta \le 1,\underset{\eta }{\max }\left\{\frac{20+\eta }{12}\right\}=1.75,$, which completes the proof.     □

We will use the following error bound lemma for spline interpolation by Hall [39]. For π : a = x₀ < x₁ <, …, < x_N = b, let $\underset{\_}{h}=\underset{j}{\min }\{x_j-x_{j-1}\},\overline{h}=\underset{j}{\max }\{x_j-x_{j-1}\}$, and $\nu =\overline{h}/\underset{\_}{h}$. For simplicity, we write $U^{i+1}\left(x\right)=\tilde{U}\left(x\right)$ and ${\Phi }^i\left(x\right)=\tilde{\Phi }$

Lemma 4.2 ([39]). Let Y be the cubic spline associate with $\tilde{U}\in C^4\left[a,b\right]$ and the partitioning π. Then,

$\begin{array}{l}\Vert {\tilde{U}}^{\left(r\right)}-Y^{\left(r\right)}\Vert \le {\lambda }_r\Vert {\tilde{U}}^{\left(4\right)}\Vert {\left(\overline{h}\right)}^{4-r},r=0,1,2,\end{array}$

where ${\lambda }_0=5/384,{\lambda }_1=\left(1/216\right)\left(9+\sqrt{3}\right),{\lambda }_2=\left(1/12\right)\left(3\nu +1\right)$

Since the mesh we have used is uniform, $\underset{\_}{h}=\overline{h}=h$.

Theorem 4.1. Let S(x, η) be the collocation approximation from the space of cubic spline ${\hat{Q}}_3\left({\Omega }_x^N\right)$ to the solution $\tilde{U}\left(x\right)$ of ordinary differential equation. If $\tilde{\Phi }\left(x_j\right)\in C^2\left[0,1\right]$, the parameter uniform error is given by

$\begin{array}{l}{\displaystyle \underset{0<\epsilon \le 1}{sup}}{\displaystyle \underset{0\le j\le N}{max}}|\hat{U}\left(x_j\right)-S\left(x_j,\eta \right)|\le C\frac{h^2}{\epsilon +h}\end{array}$

where h is sufficiently small and C is a positive constant independent of ε and N.

Proof. Let Y_s(x) be unique spline interpolation from ${\hat{Q}}_3\left({\Omega }_x^N\right)$ to the solution $\hat{U}\left(x\right)$ for boundary value problem Eq. (14) given by

$\begin{array}{ll}{Y}_s\left(x\right)={\displaystyle \sum _{j=-1}^{N+1}}{\overline{\delta }}_j{Q}_j\left(x\right).& (29)\end{array}$

If $\tilde{\Phi }\left(x\right)\in C^2\left(\overline{{\Omega }_x}\right)$, then $\tilde{U}\left(x\right)\in C^4\left(\overline{{\Omega }_x}\right)$, and so using Lemma 4.2, we have

$\begin{array}{ll}\Vert {\tilde{U}}^{\left(r\right)}\left(x\right)-Y_s^{\left(r\right)}\left(x\right)\Vert \le {\lambda }_r\Vert {\tilde{U}}^{\left(4\right)}\left(x\right)\Vert h^{4-r},r=0,1,2& (30)\end{array}$

where λ_r's are independent of h and N. From the estimate Eq. (30), we obtain

$\begin{array}{l}|{𝔏}_{\sigma }^{M}S\left(x_j\right)-{𝔏}_{\sigma }^M{Y}_s\left(x_j\right)|=|{𝔏}_{\epsilon }^M\hat{U}\left(x_j\right)-{𝔏}_{\sigma }^M{Y}_s\left(x_j\right)|\\ \le |{\sigma }_j-\epsilon \Vert {\tilde{U}}^{\left(2\right)}\left(x\right)|\\ +\left(\left|\epsilon \right|{\lambda }_2{h}^2+\parallel a\parallel {\lambda }_1{h}^3+\parallel d\parallel {\lambda }_0{h}^4\right)\left|{\tilde{U}}^{\left(4\right)}\left(x\right)\right|.\end{array}$

Now, from Lemma (3.4), estimate Eq. (25) and using the argument that since ε≪1 and ${\epsilon }^{-n}exp\left(-\frac{\alpha \left(1-x\right)}{\epsilon }\right)\to 0$ as ε → 0,∀x ∈ Ω_x, we easily obtain

$\begin{array}{ll}|{𝔏}_{\sigma }^{M}S\left(x_j\right)-{𝔏}_{\sigma }^M{Y}_s\left(x_j\right)|\le C\frac{h^2}{h+\epsilon }& (31)\end{array}$

Let ${𝔏}_{\sigma }^M{Y}_s\left(x_j\right)=\overline{\Phi }\left(x_j\right),\forall 0\le j\le N$ with boundary conditions Y_s(0) = ψ₀(t_i+1), Y_s(1) = ψ₁(t_i+1) leads to $R\overline{\delta }=\overline{\Phi }$. Then, it follows with Eq. (28) that

$\begin{array}{ll}R\left(\delta -\overline{\delta }\right)=\Phi -\overline{\Phi },& (32)\end{array}$

where

$\begin{array}{l}\delta -\overline{\delta }={\left({\delta }_0-{\overline{\delta }}_0,{\delta }_1-{\overline{\delta }}_1,\dots ,{\delta }_N-{\overline{\delta }}_N\right)}^t\\ \Phi -\overline{\Phi }={\left(h^2\left({\Phi }_0^i-{\overline{\Phi }}_0^i\right),h^2\left({\Phi }_1^i-{\overline{\Phi }}_1^i\right),\dots ,h^2\left({\Phi }_N^i-{\overline{\Phi }}_N^i\right)\right)}^t\hfill \end{array}$

It is obvious from inequality Eq. (31) that

$\begin{array}{ll}\Phi -\overline{\Phi }=C\frac{h^4}{h+\epsilon }& (33)\end{array}$

In Eq. 19, $\left|R_j^0\right|>\left\{\left|R_j^{-}\right|+\left|R_j^{+}\right|\right\}$ implies $\left(2+\eta \right)\left(2{\sigma }_j\pm h^2\frac{d_j^{i+1}}{6}\right)>0$. Therefore for η > −2 and sufficiently small values of h, matrix R is strictly diagonally dominant, thus non singular. So by the estimate given [40], we get

$\begin{array}{ll}\parallel R^{-1}\parallel \le \frac{C}{h^2}& (34)\end{array}$

Combining of Eq. (32)–Eq. (34) gives

$\begin{array}{ll}|\delta -\overline{\delta }|\le C\frac{h^2}{h+\epsilon }& (35)\end{array}$

Using Eq. (26) and Eq. (27), we have

$|{\delta }_{-1}-{\overline{\delta }}_{-1}|\le C\frac{h^2}{h+\epsilon }\text{and }|{\delta }_{N+1}-{\overline{\delta }}_{N+1}|\le C\frac{h^2}{h+\epsilon }$

Therefore,

$\begin{array}{ll}{\displaystyle \underset{-1\le j\le N+1}{max}}|{\delta }_j-{\overline{\delta }}_j|\le C\frac{h^2}{h+\epsilon }& (36)\end{array}$

We define $e_j^{*}=|{\delta }_j-{\overline{\delta }}_j|$. Then,

$e^{*}={\displaystyle \underset{-1\le j\le N+1}{max}}\left\{e_j^{*}\right\}\to 0\text{as}h\to 0.$

This proves that the proposed finite difference scheme is unconditionally stable. The above inequality Eq. (36) together with Lemma 4.1 enables us to estimate |S(x)−Y_s(x)|; hence, $|Y_s\left(x\right)-\tilde{U}\left(x\right)|$. We have

$\begin{array}{l}|S\left(x\right)-Y_s\left(x\right)|=\left|{\displaystyle \sum _{j=-1}^{N+1}}\left({\delta }_j-{\overline{\delta }}_j\right)Q_j\left(x,\eta \right)\right|\\ \le {\displaystyle \underset{-1\le j\le N+1}{max}}|{\delta }_j-{\overline{\delta }}_j|{\displaystyle \sum _{j=-1}^{N+1}}\left|Q_j\left(x,\eta \right)\right|\\ \le C\frac{h^2}{h+\epsilon }\end{array}$

This gives $\underset{0\le j\le N}{max}|S\left(x_j\right)-Y_s\left(x_j\right)|\le C\frac{h^2}{h+\epsilon }$. Therefore, using the triangular inequality, we obtain

${\displaystyle \underset{0<\epsilon \le 1}{sup}}{\displaystyle \underset{0\le j\le N}{max}}|\hat{U}\left(x_j\right)-S\left(x_j,\eta \right)|\le C\frac{h^2}{\epsilon +h},$

hence the result.     □

Theorem 4.2. Let u(x, t) be the solution of problem Eq. (1) and S(x, t_i+1) be the collocation approximation from the space ${\hat{Q}}_3\left({\Omega }_x^N\right)$ to the solution u(x, t_i+1) at the (i + 1)^th time level of the fully discretized scheme after the temporal discretization. If $\Phi \left(x,t_{i+1}\right)\in C^2\left({\Omega }_x\right)$, the uniform error estimate is given by

$\begin{array}{l}\parallel u\left(x_j,t_{i+1}\right)-S\left(x_j,t_{i+1}\right)\parallel \le C\left(\Delta t+\frac{h^2}{\epsilon +h}\right)\end{array}$

Proof. The proof is the consequence of Lemma 3.3 and Theorem 4.1.     □

5. Numerical examples, results, and discussions

In this section, we will see the applicability of the proposed method by considering two test problems. For each η = −1 value in the range [−8, 1], the computation was carried out to find the most valuable free parameter, which will give a minimum error. The least absolute errors were found for η = −1. For the remainder of the calculation, we have set up η = −1.

Example 5.1 ([41]). Consider the following test problem

$\begin{array}{c}\frac{\partial u\left(x,t\right)}{\partial t}-\epsilon \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+\left(2-x^2\right)\frac{\partial u\left(x,t\right)}{\partial x}+xu\left(x,t\right)\hfill \\ =-u\left(x,t-1\right)+10t^{2}exp\left(-t\right)x\left(1-x\right),\left(x,t\right)\in {\Omega }_x\times {\Omega }_t,\hfill \\ u\left(x,t\right)=0,\left(x,t\right)\in \overline{{\Omega }_x}\times {\Omega }_{t*},\hfill \\ u\left(0,t\right)=0,u\left(1,t\right)=0,t\in \overline{{\Omega }_t}\hfill \end{array}$

We employ the following double mesh principle to determine the absolute error in the solution because the test problem's exact solution is unknown. For each ε, the maximum point wise error is given as

$\begin{array}{ll}{E}_{\epsilon }^{N,M}={\displaystyle \underset{i}{max}}\left({\displaystyle \underset{j}{max}}|S_{x_{2j},t_{2i}}^{2N,2M}-S_{x_j,t_i}^{N,M}|\right)& (37)\end{array}$

where S^{N, M} and S^{2N, 2M} are the computed solutions obtained on two different meshes D^{N, M} and D^{2N, 2M} respectively. D^{2N, 2M} obtained from D^{N, M} by the interpolation technique. The corresponding order of convergence is computed as

$\begin{array}{ll}{P}_{\epsilon }^{N,M}={log}_2\left(\frac{E_{\epsilon }^{N,M}}{E_{\epsilon }^{2N,2M}}\right)& (38)\end{array}$

The ε-uniform point wise error E^{N, M} is estimated as

$\begin{array}{ll}{E}^{N,M}={\displaystyle \underset{\epsilon }{max}}{E}_{\epsilon }^{N,M}& (39)\end{array}$

and the ε-uniform order of convergence is calculated as

$\begin{array}{ll}{P}^{N,M}={log}_2\left(\frac{E^{N,M}}{E^{2N,2M}}\right)& (40)\end{array}$

Example 5.2. Consider the following problem

$\begin{array}{c}\frac{\partial u\left(x,t\right)}{\partial t}-\epsilon \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+\left(1+x\left(1-x\right)\right)\frac{\partial u\left(x,t\right)}{\partial x}\hfill \\ =-u\left(x,t-1\right)+f\left(x,t\right),\left(x,t\right)\in {\Omega }_x\times {\Omega }_t,\hfill \\ u\left(x,t\right)={\psi }_b\left(x,t\right),\left(x,t\right)\in \overline{{\Omega }_x}\times {\Omega }_{t*},\hfill \\ u\left(0,t\right)=0,u\left(1,t\right)=0,t\in \overline{{\Omega }_t}\hfill \end{array}$

We choose the initial data ψ_b(x, t) and the source function f(x, t) to fit with the exact solution

$\begin{array}{l}u\left(x,t\right)=exp\left(-t\right)\left\{p_1+p_{2}x-exp\left(-\left(1-x\right)/\epsilon \right)\right\},\hfill \end{array}$

where p₁ = exp(−1/ε) and p₂ = 1−p₁.

As we know the exact solution, we compute the point-wise error as

$\begin{array}{ll}{\tilde{E}}_{\epsilon }^{N,M}={\displaystyle \underset{i}{max}}\left({\displaystyle \underset{j}{max}}|u_{x_j,t_i}^{N,M}-S_{x_j,t_i}^{N,M}|\right),& (41)\end{array}$

where $u_{x_j,t_i}^{N,M}$ and $S_{x_j,t_i}^{N,M}$ denote the exact and numerical solution obtained on D^{N, M}. The corresponding computed order of convergence is calculated as

$\begin{array}{ll}{\tilde{P}}_{\epsilon }^{N,M}={log}_2\left(\frac{{\tilde{E}}_{\epsilon }^{N,M}}{{\tilde{E}}_{\epsilon }^{2N,2M}}\right)& (42)\end{array}$

The ε- uniform point wise error ${\tilde{E}}^{N,M}$ and the corresponding order of convergence ${\tilde{P}}^{N,M}$ are calculated as

$\begin{array}{ll}{\tilde{E}}^{N,M}={\displaystyle \underset{\epsilon }{max}}{\tilde{E}}_{\epsilon }^{N,M}& (43)\end{array}$

$\begin{array}{ll}{\tilde{P}}^{N,M}={log}_2\left(\frac{{\tilde{E}}^{N,M}}{{\tilde{E}}^{2N,2M}}\right)& (44)\end{array}$

The numerical solution profiles for Example 5.1 and Example 5.2 at N = 64 and different values ε are graphically represented in Figures 1, 2, respectively. These figures shows the existence of boundary layer at x = 1, and it is clearly observed that width of the boundary layer decreases as perturbation parameter ε decreases, which is the effect of ε. Figures 3, 4 provide the numerical solution for test problem Example 5.1 and Example 5.2 for different values of t for fixed value of perturbation parameter ε at η = −1, where N = 64. Figure 5 provides the log-log plot of maximum absolute errors for Example 5.1 and Example 5.2. Figure 6 represents the graphs of the solution for various values of ε. Since the plots follow a straight line, this shows us that the maximum absolute point-wise error changes as a constant power of mesh parameter N. In addition, the negative slope of the lines states that the maximum absolute error decreases as the number of mesh points increases. In these figures, the plots are parallel, which shows the parameter-uniform convergence of the scheme. The maximum point wise error $(E_{\epsilon }^{N,M},{\tilde{E}}_{\epsilon }^{N,M})$, ε-uniform errors $(E^{N,M},{\tilde{E}}^{N,M})$, rate of convergence $(P_{\epsilon }^{N,M},{\tilde{P}}_{\epsilon }^{N,M})$, and ε-uniform rate of convergence $(P^{N,M},{\tilde{P}}^{N,M})$ for Example 5.1 and Example 5.2 are presented in Tables 2, 3, respectively, at η = −1. The numerical results presented in Tables 2, 3 show that the proposed method is ε-uniformly convergent as for fixed value of ε. When N and M are increases, the maximum point wise error $(E_{\epsilon }^{N,M},{\tilde{E}}_{\epsilon }^{N,M})$ and the maximum nodal errors $(E^{N,M},{\tilde{E}}^{N,M})$ decreases. We see that the maximum point wise error and the rate of convergence stabilize as ε → 0 for each N and M. A comparison of maximum point wise error for 5.1 calculated by the proposed method at η = −1 and for in Kumar [8] is presented in Table 4. Computational results in Table 4 shows the proposed method provide more accurate solutions than in Kumar [8]. Furthermore, we note that all computations have been performed using MATLAB^® R2022b software package (The Mathworks, Inc.), on a 64 bit Windows 11 hp CPU PC machine, with Intel(R) Core(TM) i3-3110M processor running at 2.40 GHz and 4.00 Gb RAM.

FIGURE 1

Figure 1. Numerical solution profile for Example 5.1, with N = 64, M = 80. (A) ε = 2⁻⁴. (B) ε = 2⁻¹⁵.

FIGURE 2

Figure 2. Numerical solution profile for Example 5.2, with N = M = 64. (A) ε = 2⁻⁴. (B) ε = 2⁻¹⁵.

FIGURE 3

Figure 3. Numerical solution of Example 5.1 for different values of ε and t with N = 64 and M = 80.

FIGURE 4

Figure 4. Numerical solution of Example 5.2 for different values of ε and t with N = M = 64.

FIGURE 5

Figure 5. Log-log plot of maximum point wise error (A) for Example 5.1 and (B) for Example 5.2.

FIGURE 6

Figure 6. Effect of ε on the solution behavior at T = 1.5. (A) Example 5.1. (B) Example 5.2.

TABLE 2

Table 2. Values of $E_{\epsilon }^{N,M},P_{\epsilon }^{N,M},E^{N,M},$ and P^{N, M} for Example 5.1.

TABLE 3

Table 3. ${\tilde{E}}_{\epsilon }^{N,M},{\tilde{P}}_{\epsilon }^{N,M},{\tilde{E}}^{N,M},$ and ${\tilde{P}}^{N,M}$ for Example 5.2.

TABLE 4

Table 4. Comparison of the maximum point-wise error for Example 5.1.

6. Conclusion

In this study, we provided a parameter uniform numerical scheme is developed to solve singularly perturbed parabolic convection-diffusion initial boundary value problems with large delay. The method is based on the implicit Euler method for temporal discretization and extended cubic B-spline collocation method with a blending function of degree four for spatial discretization using artificial viscosity both on the uniform mesh. The theoretical results which show the parameter-uniform convergence of the method are established and the proposed method is shown order $\mathcal{O}\left(\Delta t+\frac{h^2}{\epsilon +h}\right)$. The scheme is also unconditionally stable. The appropriate choice of the free parameter η minimizes the error. To validate the theoretical results two test examples are presented. Graphical and tabular representations of the solutions and accuracy of the examples' results are provided. The numerical results obtained by the proposed method are compared with the numerical results in some existing literature.

As future directions of this study, we extend the proposed scheme for solving non-linear and higher dimensional singularly perturbed delay partial differential equations with Dirichlet boundary conditions.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

ZH: Conceptualization, Formal analysis, Investigation, Methodology, Resources, Software, Writing—original draft, Writing—review and editing. GD: Conceptualization, Formal analysis, Investigation, Methodology, Resources, Software, Supervision, Writing—original draft, Writing—review and editing. All authors contributed to the article and approved the submitted version.

Acknowledgments

The authors kindly thank the editor for handling the manuscript carefully and referees for their valuable contributions in modifying the original version of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Cooke K, Van den Driessche P, Zou X. Interaction of maturation delay and nonlinear birth in population and epidemic models. J Math Biol. (1999) 39:332–52. doi: 10.1007/s002850050194

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Wang XT. Numerical solution of delay systems containing inverse time by hybrid functions. Appl Math Comput. (2006) 173:535–46. doi: 10.1016/j.amc.2005.04.056

CrossRef Full Text | Google Scholar

3. Zhao T. Global periodic-solutions for a differential delay system modeling a microbial population in the chemostat. J Math Anal Appl. (1995) 193:329–52. doi: 10.1006/jmaa.1995.1239

CrossRef Full Text | Google Scholar

4. Smolen P, Baxter DA, Byrne JH. A reduced model clarifies the role of feedback loops and time delays in the Drosophila circadian oscillator. Biophys J. (2002) 83:2349–59. doi: 10.1016/S0006-3495(02)75249-1

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Vielle B, Chauvet G. Delay equation analysis of human respiratory stability. Math Biosci. (1998) 152:105–22. doi: 10.1016/S0025-5564(98)10028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Villasana M, Radunskaya A. A delay differential equation model for tumor growth. J Math Biol. (2003) 47:270–94. doi: 10.1007/s00285-003-0211-0

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wang PKC. Asymptotic stability of a time-delayed diffusion system. J Appl Mech. (1963) 30:500–4. doi: 10.1115/1.3636609

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kumar D. A parameter-uniform scheme for the parabolic singularly perturbed problem with a delay in time. Numer Methods Partial Differ Equ. (2021) 37:626–42. doi: 10.1002/num.22544

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wu J. Theory and Applications of Partial Functional Differential Equations, Vol 119. New York, NY: Springer Science & Business Media (1996). doi: 10.1007/978-1-4612-4050-1

CrossRef Full Text | Google Scholar

10. Kuang Y. Delay Differential Equations: with Applications in Population Dynamics. Cambridge, MA: Academic press (1993).

Google Scholar

11. Das A, Natesan S. Uniformly convergent hybrid numerical scheme for singularly perturbed delay parabolic convection-diffusion problems on Shishkin mesh. Appl Math Comput. (2015) 271:168–86. doi: 10.1016/j.amc.2015.08.137

CrossRef Full Text | Google Scholar

12. Negero N, Duressa G. An efficient numerical approach for singularly perturbed parabolic convection-diffusion problems with large time-lag. J Math Model. (2022) 10:173–110. doi: 10.1016/j.rinam.2022.100338

CrossRef Full Text | Google Scholar

13. Negero NT, Duressa GF. A method of line with improved accuracy for singularly perturbed parabolic convection-diffusion problems with large temporal lag. Results Appl Math. (2021) 11:100174. doi: 10.1016/j.rinam.2021.100174

CrossRef Full Text | Google Scholar

14. Woldaregay MM, Aniley WT, Duressa GF. Novel numerical scheme for singularly perturbed time delay convection-diffusion equation. Adv Math Phys. (2021) 2021:1–13. doi: 10.1155/2021/6641236

CrossRef Full Text | Google Scholar

15. Kumar K, Podila PC, Das P, Ramos H. A graded mesh refinement approach for boundary layer originated singularly perturbed time-delayed parabolic convection diffusion problems. Math Methods Appl Sci. (2021) 44:12332–50. doi: 10.1002/mma.7358

CrossRef Full Text | Google Scholar

16. Abdelhakem M, Youssri Y. Two spectral Legendre's derivative algorithms for Lane-Emden, Bratu equations, and singular perturbed problems. Appl Numer Math. (2021) 169:243–55. doi: 10.1016/j.apnum.2021.07.006

CrossRef Full Text | Google Scholar

17. Abd-Elhameed W, Youssri Y, Doha E. A novel operational matrix method based on shifted Legendre polynomials for solving second-order boundary value problems involving singular, singularly perturbed and Bratu-type equations. Math Sci. (2015) 9:93–102. doi: 10.1007/s40096-015-0155-8

CrossRef Full Text | Google Scholar

18. Salama A, Al-Amery D. A higher order uniformly convergent method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2017) 94:2520–46. doi: 10.1080/00207160.2017.1284317

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

20. Babu G, Bansal K. A high order robust numerical scheme for singularly perturbed delay parabolic convection diffusion problems. J Appl Math Comput. (2022) 68:363–89. doi: 10.1007/s12190-021-01512-1

CrossRef Full Text | Google Scholar

21. Govindarao L, Mohapatra J. A second order numerical method for singularly perturbed delay parabolic partial differential equation. Eng Comput. (2018) 36:420–44. doi: 10.1108/EC-08-2018-0337

CrossRef Full Text | Google Scholar

22. Podila PC, Kumar K. A new stable finite difference scheme and its convergence for time-delayed singularly perturbed parabolic PDEs. Comput Appl Math. (2020) 39:1–16. doi: 10.1007/s40314-020-01170-2

CrossRef Full Text | Google Scholar

23. Hailu WS, Duressa GF. Parameter-uniform cubic spline method for singularly perturbed parabolic differential equation with large negative shift and integral boundary condition. Res Math. (2022) 9:2151080. doi: 10.1080/27684830.2022.2151080

CrossRef Full Text | Google Scholar

24. Hassen ZI, Duressa GF. New approach of convergent numerical method for singularly perturbed delay parabolic convection-diffusion problems. Res Math. (2023) 10:2225267. doi: 10.1080/27684830.2023.2225267

CrossRef Full Text | Google Scholar

25. Daba IT, Duressa GF. Extended cubic B-spline collocation method for singularly perturbed parabolic differential-difference equation arising in computational neuroscience. Int J Numer Method Biomed Eng. (2021) 37:e3418. doi: 10.1002/cnm.3418

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Daba IT, Duressa GF. Collocation method using artificial viscosity for time dependent singularly perturbed differential-difference equations. Math Comput Simul. (2022) 192:201–20. doi: 10.1016/j.matcom.2021.09.005

CrossRef Full Text | Google Scholar

27. Kumar D, Kadalbajoo MK. A parameter-uniform numerical method for time-dependent singularly perturbed differential-difference equations. Appl Math Model. (2011) 35:2805–19. doi: 10.1016/j.apm.2010.11.074

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

29. Kumar D, Kumari P. A parameter-uniform scheme for singularly perturbed partial differential equations with a time lag. Numer Methods Partial Differ Equ. (2020) 36:868–86. doi: 10.1002/num.22455

CrossRef Full Text | Google Scholar

30. Stynes M, Stynes D. Convection-diffusion Problems: An Introduction to their Analysis and Numerical Solution, Vol 196. Washington, DC: American Mathematical Society; Atlantic Association for Research in the Mathematical Sciences (2018). doi: 10.1090/gsm/196

CrossRef Full Text | Google Scholar

31. Mbroh NA, Noutchie SCO, Massoukou RYM. A robust method of lines solution for singularly perturbed delay parabolic problem. Alex Eng J. (2020) 59:2543–54. doi: 10.1016/j.aej.2020.03.042

CrossRef Full Text | Google Scholar

32. Rothe E. Zweidimensionale parabolische randwertaufgaben als grenzfall eindimensionaler randwertaufgaben. Math Ann. (1930) 102:650–70. doi: 10.1007/BF01782368

CrossRef Full Text | Google Scholar

33. Kadalbajoo MK, Gupta V, Awasthi A. A uniformly convergent B-spline collocation method on a nonuniform mesh for singularly perturbed one-dimensional time-dependent linear convection-diffusion problem. J Comput Appl Math. (2008) 220:271–89. doi: 10.1016/j.cam.2007.08.016

CrossRef Full Text | Google Scholar

34. Kellogg RB, Tsan A. Analysis of some difference approximations for a singular perturbation problem without turning points. Math Comput. (1978) 32:1025–39. doi: 10.1090/S0025-5718-1978-0483484-9

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Clavero C, Jorge J, Lisbona F. A uniformly convergent scheme on a nonuniform mesh for convection-diffusion parabolic problems. J Comput Appl Math. (2003) 154:415–29. doi: 10.1016/S0377-0427(02)00861-0

CrossRef Full Text | Google Scholar

36. Dağ İ, Irk D, Sarı M. The extended cubic B-spline algorithm for a modified regularized long wave equation. Chin Phys B. (2013) 22:040207. doi: 10.1088/1674-1056/22/4/040207

CrossRef Full Text | Google Scholar

37. Gang X, Guo-Zhao W. Extended cubic uniform B-spline and α-B-spline. Acta Autom Sin. (2008) 34:980–4. doi: 10.1016/S1874-1029(08)60047-6

PubMed Abstract | CrossRef Full Text | Google Scholar

38. O'Malley RE. Singular Perturbation Methods for Ordinary Differential Equations, Vol 89. New York, NY: Springer (1991). doi: 10.1007/978-1-4612-0977-5

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Hall C. On error bounds for spline interpolation. J Approx Theory. (1968) 1:209–18. doi: 10.1016/0021-9045(68)90025-7

CrossRef Full Text | Google Scholar

40. Varah JM. A lower bound for the smallest singular value of a matrix. Linear Algebra Appl. (1975) 11:3–5. doi: 10.1016/0024-3795(75)90112-3

CrossRef Full Text | Google Scholar

41. Kumar D, Kumari P. 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. doi: 10.1007/s12190-018-1174-z

CrossRef Full Text | Google Scholar