Gaussian quadrature method with exponential fitting factor for two-parameter singularly perturbed parabolic problem

Abstract

The parabolic convection-diffusion-reaction problem is examined in this work, where the diffusion and convection terms are multiplied by two small parameters, respectively. The proposed approach is based on a fitted operator finite difference method. The Crank-Nicolson method on uniform mesh is utilized to discretize the time variables in the first step. Two-point Gaussian quadrature rule is used for further discretizing these semi-discrete problems in space, and the second order interpolation of the first derivatives is utilized. The fitting factor’s value, which accounts for abrupt changes in the solution, is calculated using the theory of singular perturbations. The developed scheme is demonstrated to be second-order accurate and uniformly convergent. The proposed method’s applicability is validated by two examples, which yielded more accurate results than some other methods found in the literatures.

Introduction

This paper examines a two-parameter singularly perturbed parabolic problems(TP-SPPPs) on the domain ${\mathbb {G}}={\mathcal {Q}} \times {\mathcal {L}}$, where ${\mathcal {Q}}=(0,1)$ and ${\mathcal {L}}=(0,T]$. The boundaries of the domain is $\bigwedge =\bigwedge _{b}\cup \bigwedge _{l}\cup \bigwedge _{r}$ where $\bigwedge _{l}=\left\{ (x,t):x=0, t \in \overline{{\mathcal {L}}}\right\} , \bigwedge _{b}=\left\{ (x,t):t=0, x \in \overline{{\mathcal {Q}}} \right\}$, and $\bigwedge _{r}=\left\{ (x,t):x=1,t\in \overline{{\mathcal {L}}}\right\}$ are the left, bottom and right sides of ${\mathbb {G}}$. We consider:

$$\begin{aligned} {\left\{ \begin{array}{ll} {\mathbb {L}}_{\varepsilon _{1},\varepsilon _{2}}w\equiv & \dfrac{\partial w}{\partial t}(x,t)-\varepsilon _{1} \dfrac{\partial ^{2}w}{\partial x^{2}}(x,t)-\varepsilon _{2}a(x,t)\dfrac{\partial w}{\partial x}(x,t)+b(x,t)w(x,t)=\\ & \hspace{2cm}\varpi (x,t),\quad (x,t)\in {\mathbb {G}},\\ & w\arrowvert _{\bigwedge _{b}}=\zeta _{b}(x),\quad x \in \overline{{\mathcal {Q}}}\\ & w\arrowvert _{\bigwedge _{l}}=\zeta _{l}(t),\quad t \in \overline{{\mathcal {L}}}\\ & w\arrowvert _{\bigwedge _{r}}=\zeta _{r}(t),\quad t \in \overline{{\mathcal {L}}} \end{array}\right. } \end{aligned}$$

(1)

where $0< \varepsilon _{1}<<1$ and $0 \le \varepsilon _{2} \le 1$ are the perturbation parameters. To establish the existence and uniqueness of the solution of Eq. (1), we make the following assumptions:

• We assume that all the functions in Eq. (1) are sufficiently smooth, bounded and possess twice continuously differentiable.

• $0<\alpha<a(x,t), 0<\beta \le b(x,t)$ and $0<\varsigma \le c(x,t),\quad (x,t)\in \overline{{\mathbb {G}}}$.

• The compatibility conditions are satisfied at the corner points of the domain.

Singularly perturbed problems(SPPs) occur in numerous engineering, biological, and physical phenomena, where they are often described by boundary layer problems involving various types of partial differential equations(PDEs) [6]. SPPs with two small parameters are used to model phenomena such as chemical reactor theory [4], flow through unsaturated porous media [5], lubrication theory [6, 7], 1D vertical groundwater hydrology and transport phenomena appearing in chemistry and biology [8].

The fact that many real-world applications involve SPPs emphasizes the need for this research even more. Solving fluid flow problems in channels with small aspect ratios, for instance, where the flow is characterized by the presence of boundary layers is one way to incorporate fluid dynamics into the process. This approach can be used in chemical engineering to address issues with chemical diffusion in porous media, where boundary layers are a defining feature of the diffusion process. This technique can be applied in electrical engineering to solve problems pertaining to electromagnetic wave propagation in waveguides, where boundary layers are present and contribute to the wave’s propagation [29].

Depending on the value of $\varepsilon _{2}$, Eq. (1) has different natures: (i) when $\varepsilon _{2}=1$, Eq. (1) belongs to the class of one-parameter convection-diffusion SPPs, and a boundary layer with width $O(\varepsilon )$ arises near $x=0$; (ii) when $\varepsilon _{2}=0$, Eq. (1) comes under the category of reaction-diffusion SPPs, and two boundary layers, each with width $O(\sqrt{\varepsilon })$, are observed at both ends $x=0$ and $x=1$. (iii) With the exception of (i) and (ii), that is, $\varepsilon _{2} \ne 0,1$, we have TP-SPPPs for which two distinct cases, $\varepsilon _{2}^{2}/\varepsilon _{1}<<1$ and $\varepsilon _{2}^{2}/\varepsilon _{2}>>1$, are contributed by the ratio of $\varepsilon _{2}^{2}$ to $\varepsilon _{1}$ and two boundary layers are observed at $x=0$ and $x=1$(see [9,10,11]).

The TP-SPPs have been treated numerically by many authors, in particular, several numerical studies have been made for TP-SPPs(see [1,2,3, 16, 18]). In the context of PDEs, however, numerous numerical studies of TP-SPPs have been conducted. To be specific, the problem under study in this work has not been thoroughly explored, even though some numerical approaches for solving these problems already exist. For instance, Jha and Kadalbajoo [12] used the implicit Euler method for time and the upwind scheme on Shishkin mesh for spatial discretization to develop a finite difference method(FDM) for TP-SPPPs. Shivhare et al. [13], developed second order accurate uniformly convergent schemes using quadratic B-splines on exponentially graded mesh. Zahra et al. [14], constructed a numerical method comprising the Euler scheme with an equidistant mesh in the temporal discretization and the exponential spline scheme with a predefined non-uniform mesh in the spatial direction for the solution of 1D SPPs. To find more literature on this, the redear may refer( [9,10,11, 15, 17, 19,20,21]).

The aforementioned studies serve as our inspiration for proposing and analyzing the parameter-uniform numerical solution for TP-SPPs. To the best of the author’s understanding, no existing numerical method in the literature incorporates two-point Gaussian quadrature rule for TP-SPPs. An exponentially fitted operator strategy for solving a SPPs with a right boundary layer is presented in [22]. In this work, we first apply the Crank-Nicolson method on a uniform mesh to discretize the time variables. Next, we employ a Taylor series expansion to transform the second-order differential equation(DE) into a first-order delay DE with a small deviating argument $\varepsilon _{1}$. Finally, we utilize the two-point Gaussian quadrature rule and the FDM to approximate the first-order DE, resulting in a fitting factor and a tridiagonal system of equations.

The Gaussian quadrature rule operates by utilizing a weighted sum of function values at particular points, referred to as nodes, to approximate the integral of the function. The purpose of selecting these nodes and weights was to reduce the approximation’s error. When dealing with problems that are singularly perturbed, the nodes and weights are selected to account for boundary layers and to offer high accuracy in the area where the solution varies quickly. When it comes to solving SPPs, the two-point Gaussian quadrature rule is an effective tool that provides superb accuracy and stability even when other numerical techniques are unable to do so.

The rest of the paper is organized as follows: We examine properties of the continuous problem in section . The semi-discretized and fully discretized scheme of the continuous problem is explained in section . In section , the parameter-uniform convergence analysis is presented. Section  presents the numerical results, and Section  summarizes the main conclusions of the paper.

Solution bounds for continuous problem

For the purpose of conducting a convergence analysis of the suggested approach, we will derive derivative bounds for the solution using the maximum principle for ${\mathbb {L}}_{\varepsilon _{1},\epsilon _{2}}w$.

Lemma 1

(Maximum Principle). Let $\varXi (x,t)\in C^{2}({\mathbb {G}})\cap C^{0}(\overline{{\mathbb {G}}})$, such that $\varXi (x,t)\ge 0, \forall (x,t)\in \bigwedge =\bigwedge _{b}\cup \bigwedge _{l}\cup \bigwedge _{r}$ & ${\mathbb {L}}_{\varepsilon _{1},\varepsilon _{2}}\varXi (x,t)\ge 0$ then, $\varXi (x,t)\ge 0, \forall (x,t)\in \overline{{\mathbb {G}}}$.

Proof

The proof is given in [13, 20]. $\square$

Lemma 2

(Uniform Stability Estimate). The solution w(x, t) of Eq. (1) is bounded and satisfies the estimate below

$$\begin{aligned} \left\| w(x,t) \right\| \le \left| \xi ^{-1}\right| \left\| \varpi (x,t) \right\| +\max \left\{ \left\| \zeta _{b}(x)\right\| ,\left\| \zeta _{l}(t)\right\| ,\left\| \zeta _{r}(t)\right\| \right\} , \forall (x,t) \in \overline{{\mathbb {G}}}, \end{aligned}$$

(2)

where $0<\xi \le b(x,t)$.

Proof

The proof is given in [13, 20]. $\square$

Lemma 3

For $\iota ,m \in \mathbb {Z^{+}}$ satisfying $0 \le \iota +2m\le 4$, the solution w(x, t) and its derivatives of Eq. (1), satisfy the bound:

$$\begin{aligned} \left\| \dfrac{\partial ^{\iota +m}w}{\partial x^{\iota }\partial t^{m}}\right\| \le \left\{ \begin{array}{rl} C\varepsilon _{1}^{-\frac{\iota }{2}}& \hbox { for}\ \varepsilon _{2}^{2}\le \frac{\varepsilon _{1}\lambda }{\alpha }, \\ C\left( \frac{\varepsilon _{2}}{\varepsilon _{1}} \right) ^{\iota }\left( \frac{\varepsilon _{2}^{2}}{\varepsilon _{1}} \right) ^{m} & \hbox { for}\ \varepsilon _{2}^{2}\ge \frac{\varepsilon _{1}\lambda }{\alpha }, \end{array} \right. \end{aligned}$$

(3)

where $\lambda \approx \min _{(x,t)\in {\overline{\varDelta }}}\dfrac{b(x,t)}{a(x,t)}$ and C is independent of $\varepsilon _{1}$ and $\varepsilon _{2}$.

Proof

For the proof of Lemma (3), refer [13]. $\square$

Construction of numerical methods

Temporal semi-discretization

The continuous problem Eq. (1) is semi-discretized in the time direction on the domain of time interval [0, T] using Crank-Nicolson method with uniform mesh. The equidistant mesh is generated using a time mesh step size of $k=T/M$

$$\begin{aligned} \Lambda ^{M}_{t}=\left\{ t_{j}=kj, j=0,1,...,M, \quad t_{M}=T, k=\frac{T}{M}\right\} , \end{aligned}$$

where M is the number of mesh elements in the time direction. The Tayloy series expansion was utilized for w(x, t) at the point $\left( x,t_{j+1/2} \right)$ to derive the Crank-Nicolson scheme in temporal direction, resulting in:

$$\begin{aligned} w^{j+1}(x)= & w^{j+1/2}(x)+\frac{k}{2}\dfrac{dw^{j+1/2}(x)}{dt}+ \frac{k^{2}}{8}\dfrac{d^{2}w^{j+1/2}(x)}{dt^{2}}+ \frac{k^{3}}{48}\dfrac{d^{3}w^{j+1/2}(x)}{dt^{3}}+... \end{aligned}$$

(4)

$$\begin{aligned} w^{j}(x)= & w^{j+1/2}(x)-\frac{k}{2}\dfrac{dw^{j+1/2}(x)}{dt}+ \frac{k^{2}}{8}\dfrac{d^{2}w^{j+1/2}(x)}{dt^{2}}- \frac{k^{3}}{48}\dfrac{d^{3}w^{j+1/2}(x)}{dt^{3}}+... \end{aligned}$$

(5)

Subtracting Eq. (4) from Eq. (5), we eliminate the term $w^{j+1/2}(x)$, resulting in:

$$\begin{aligned} \dfrac{dw^{j+1/2}(x)}{dt}=\dfrac{w^{j+1}(x)-w^{j}(x)}{k}+TE^{j+1/2}(x). \end{aligned}$$

(6)

Here, $TE^{j+1/2}(x)$ is local truncation error(LTE) of the scheme and is defined as:

$$\begin{aligned} TE^{j+1/2}(x)=\frac{k^{3}}{24}\dfrac{d^{3}w^{j+1/2}(x)}{dt^{3}}+\text{ Higher } \text{ order } \text{ terms } \end{aligned}$$

(7)

which is order three. Below is the semi-discretization of Eq. (1):

$$\begin{gathered} - \varepsilon _{1} \left( {w_{{xx}} (x)} \right)^{{j + 1}} - \varepsilon _{2} a^{{j + 1}} (x)\left( {w_{x} (x)} \right)^{{j + 1}} + \left( {b^{{j + 1}} (x) + \frac{2}{k}} \right)w^{{j + 1}} (x)\, = \, \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\varepsilon _{1} \left( {w_{{xx}} (x)} \right)^{j} + \varepsilon _{2} a^{j} (x)\left( {w_{x} (x)} \right)^{j} - \left( {b^{j} (x) - \frac{2}{k}} \right)w^{j} (x) + \varpi ^{{j + 1}} (x) + \varpi ^{j} (x). \hfill \\ \end{gathered}$$

(8)

At the $(j+1)th$ time step, the semi-discrete solution $w(x,t_{j+1})$ of Eq. (1) satisfies the following differential equation of the space variable x.

$$\left\{ {\begin{array}{*{20}l} {\mathbb{L}_{{\varepsilon _{1} ,\varepsilon _{2} }}^{M} w^{{j + 1}} (x) = - \varepsilon _{1} \left( {w_{{xx}} (x)} \right)^{{j + 1}} - \varepsilon _{2} a^{{j + 1}} (x)\left( {w_{x} (x)} \right)^{{j + 1}} + } \hfill \\ {\left( {b^{{j + 1}} (x) + \frac{2}{k}} \right)w^{{j + 1}} (x) = {\mathcal{H}}^{j} (x),} \hfill \\ {w^{{j + 1}} (0) = \zeta _{l} (t_{{j + 1}} ),\quad w^{{j + 1}} (1) = \zeta _{r} (t_{{j + 1}} ),\quad 0 \le j \le M,} \hfill \\ {w^{0} (x) = \zeta _{b} (x),\quad x \in (0,1),} \hfill \\ \end{array} } \right.$$

(9)

where

$$\begin{aligned} {\mathcal {H}}^{j}(x)=\varepsilon _{1}\left( w_{xx}(x)\right) ^{j}+\varepsilon _{2} a^{j}(x)\left( w_{x}(x)\right) ^{j}-\left( b^{j}(x) -\dfrac{2}{k}\right) w^{j}(x)+\varpi ^{j+1}(x)+\varpi ^{j}(x). \end{aligned}$$

Lemma 4

(Semi-discrete Maximum Principle.) At $(j+1)$th time level, let $\varPsi (x,t_{j+1})\in C^{2}(\overline{{\mathcal {Q}}})$ be smooth function such that $\varPsi (0,t_{j+1})\ge 0$ and $\varPsi (1,t_{j+1})\ge 0$, then ${\mathbb {L}}_{\varepsilon _{1},\varepsilon _{2}}^{M}\varPsi (x,t_{j+1})\ge 0, \forall x\in {\mathbb {G}}$ implies $\varPsi (x,t_{j+1})\ge 0, \forall x \in \overline{{\mathcal {Q}}}$.

Proof

Let $\left( \chi ^{*},t_{j+1}\right) \in \left\{ \left( x,t_{j+1}\right) :x\in \overline{{\mathcal {Q}}}\right\}$ be such that $\varPsi (\chi ^{*},t_{j+1})=\min _{x\in \overline{{\mathcal {Q}}}}\varPsi (x,t_{j+1})$, and let $\varPsi (\chi ^{*},t_{j+1})<0$. Then $\left( \chi ^{*},t_{j+1}\right) \notin \left\{ (0,t_{j+1}),(1,t_{j+1}) \right\}$, also $\varPsi _{x}(\chi ^{*},t_{j+1})=0$ and $\varPsi _{xx}(\chi ^{*},t_{j+1})\ge 0$. By using this assumption and Eq. (9), we arrive ${\mathbb {L}}_{\varepsilon _{1},\varepsilon _{2}}^{M}\varPsi (x,t_{j+1})<0$, which is a contradiction to the hypothesis, ${\mathbb {L}}_{\varepsilon ,\mu }^{M}\varPsi (x,t_{j+1})\ge 0$. Hence, our asumption is wrong so $\varPsi (\chi ^{*},t_{j+1})\ge 0$ & this implies that $\varPsi (x,t_{j+1})\ge 0, \forall x \in \overline{{\mathcal {Q}}}$. $\square$

Definition 1

[23] The global and local errors of the time semi-discretization method in Eq. (9) are given by

$${\mathcal{E}}^{{j + 1}} = \max _{{j + 1 \le T/k}} \left\| {w^{{j + 1}} (x) - W^{{j + 1}} (x)} \right\|,\quad e^{{j + 1}} = \left| {w^{{j + 1}} (x) - W^{{j + 1}} (x)} \right|,$$

(10)

The local error quantifies each time steps contribution to the overall error of time semi-discretization. Instead of using $W^{j+1}(x)$ as initial data, the exact values $w^{j+1}(x)$ are used to calculate LTE after each time semi-discretization step. By applying the lemmas provided, it is possible to establish the consistency and stability of the Crank-Nicolson method, which are essential for determining the convergence rate in the temporal domain.

Lemma 5

If the solution w(x, t) of Eq. (1) and its time derivatives are bounded in $(x,t) \in [0,1]\times [0,T]$, independent of $\varepsilon _{1},\varepsilon _{2},N$ & M, then the LTE associated with Eq. (9) satisfies $\left\| e^{j+1} \right\| \le C(k^{3})$.

Proof

The proof is directily from Eq. (6) and Eq. (7). $\square$

Lemma 6

The global error estimate ${\mathcal {E}}^{j+1}$ at the $(j+1)th$ time level is bounded by ${\mathcal {E}}^{j+1}\le C(k^{2}),\forall j\le T/k$, with a constant C that is independent of $\varepsilon _{1},\varepsilon _{2},x$ & t.

Proof

The complete proof can be found in [17, 24] $\square$

To grasp the widths of the boundary layers in Eq. (9), consider the characteristic equation:

$$\begin{aligned} \varepsilon _{1} \lambda ^{2}(x)+ \varepsilon _{2}a^{j+1}(x)\lambda (x)-\beta ^{j+1}(x)=0, \end{aligned}$$

(11)

where $\beta ^{j+1}(x)=b^{j+1}(x)+2/k$. The solution of Eq. (11) establish two continuous functions:$\lambda _{1,2}(x):[0,1]\rightarrow {\mathbb {R}}$, where

$$\begin{gathered} \lambda _{1} (x) = - \frac{{\varepsilon _{2} a^{{j + 1}} (x)}}{{2\varepsilon _{1} }} + \sqrt {\left( {\frac{{\varepsilon _{2} a^{{j + 1}} (x)}}{{2\varepsilon _{1} }}} \right)^{2} + \frac{{\beta ^{{j + 1}} (x)}}{{\varepsilon _{1} }}} \ge 0, \hfill \\ \lambda _{2} (x) = - \frac{{\varepsilon _{2} a^{{j + 1}} (x)}}{{2\varepsilon _{1} }} - \sqrt {\left( {\frac{{\varepsilon _{2} a^{{j + 1}} (x)}}{{2\varepsilon _{1} }}} \right)^{2} + \frac{{\beta ^{{j + 1}} (x)}}{{\varepsilon _{1} }}} \le 0. \hfill \\ \end{gathered}$$

(12)

These roots represent the boundary layer behavior of the solution in the neighborhood of $x=0$ and $x=1$(see [13, 17]). Let $\varrho _{1}=-\max _{x\in [0,1]}\lambda _{2}(x)$ and $\varrho _{2}=\min _{x\in [0,1]}\lambda _{1}(x)$, we have three cases:(i) when $\frac{\varepsilon _{2}^{2}}{\varepsilon _{1}}\rightarrow 0$, as $\varepsilon _{1} \rightarrow 0, \varrho _{1}=\varrho _{2}=\sqrt{\frac{\beta ^{j+1}(x)}{\varepsilon _{1}}}$, where $0<\beta ^{*}<\beta ^{j+1}$ and (ii) when $\frac{\varepsilon _{1}}{\varepsilon _{2}^{2}}\rightarrow 0$, as $\varepsilon _{2} \rightarrow 0, \varrho _{1}=\varrho _{2}=\frac{\varepsilon _{2} a^{*}}{\varepsilon _{1}}$ and $\varrho _{2}=0$, where $0<a^{*}<a^{j+1}$, (iii) when $\varrho _{2}<<\varrho _{1}$, the solutions exhibit a stronger boundary layer along $x=0$ and $x=1$.

Next, we give the semi-discrete bound of the solution $w^{j+1}(x)$ of Eq. (9).

Lemma 7

The solution $w^{j+1}(x)$, of Eq. (9) satisfies the derivative bound

$$\left| {\frac{{d^{q} w^{{j + 1}} (x)}}{{dx^{q} }}} \right| \le C\left( {1 + {\varrho }_{1}^{q} e^{{ - \tau {\varrho }_{1} x}} + {\varrho }_{2}^{q} e^{{ - \tau {\varrho }_{2} (1 - x)}} } \right),q = 0(1)r,r\, > \,1\,\,and\,\,0 < \,\tau < 1.$$

(13)

Proof

For the proof, see [13, 17]. $\square$

Spatial discretization

We convert the second order DE into a first order delay DE with a minor deviating argument by using a Taylor series expansion of $w^{j+1}_{x}(x+\varepsilon _{1})$.

$$w_{x}^{{j + 1}} (x + \varepsilon _{1} ) \simeq w_{x}^{{j + 1}} (x) + \varepsilon _{1} w_{{xx}}^{{j + 1}} (x) + O(\varepsilon _{1}^{2} ).$$

(14)

From Eq. (14), we have

$$\varepsilon _{1} w_{{xx}}^{{j + 1}} (x) = w_{x}^{{j + 1}} (x + \varepsilon _{1} ) - w_{x}^{{j + 1}} (x) + O(\varepsilon _{1}^{2} ).{\text{ }}$$

(15)

We obtain an asymptotically equivalent first order neutral type DE with a minor deviating argument via implementing Eq. (15) in Eq. (8):

$$\begin{gathered} - \left( {w_{x} \left( {x + \varepsilon _{1} } \right)} \right)^{{j + 1}} + \left( {w_{x} (x)} \right)^{{j + 1}} - \varepsilon _{2} a^{{j + 1}} (x)\left( {w_{x} (x)} \right)^{{j + 1}} + \hfill \\ \left( {b^{{j + 1}} (x) + \frac{2}{k}} \right)w^{{j + 1}} (x) = \left( {w_{x} \left( {x + \varepsilon _{1} } \right)} \right)^{j} - \left( {w_{x} (x)} \right)^{j} + \varepsilon _{2} a^{j} \left( {w_{x} (x)} \right)^{j} \hfill \\ \,\,\,\,\, - \left( {b^{j} (x) - \frac{2}{k}} \right)w^{j} (x) + \varpi ^{{j + 1}} (x) + \varpi ^{j} (x). \hfill \\ \end{gathered}$$

(16)

Since $\varepsilon _{1}$ is small and the truncated term $\varepsilon _{1}^{2}/2$ is very small and approaches zero, the transition from Eq. (8) to Eq. (16) is allowed. Further detail on the validity of this transition can also be found in [26]. Subsequently, we discretized the spatial domain $[0,1]$ as follows: $\Lambda ^{N}_{x}=\left\{ x_{i}=ih, i=0(1)N, x_{0}=0,x_{N}=1,h=1/N \right\}$. When we integrate Eq. (16) from $\left[ x_{i-1},x_{i}\right]$, for $i=1(1)(N+1)$, we have:

$$\begin{gathered} - \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {w_{x} \left( {x + \varepsilon _{1} } \right)} \right)^{{j + 1}} } dx + \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {w_{x} (x)} \right)^{{j + 1}} } dx \hfill \\ \,\,\, - \varepsilon _{2} \int_{{x_{{i - 1}} }}^{{x_{i} }} {a^{{j + 1}} } (x)\left( {w_{x} (x)} \right)^{{j + 1}} dx + \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {b^{{j + 1}} (x) + \frac{2}{k}} \right)} w^{{j + 1}} (x)dx = \hfill \\ \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {w_{x} \left( {x + \varepsilon _{1} } \right)} \right)^{j} } dx - \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {w_{x} (x)} \right)^{j} } dx + \varepsilon _{2} \int_{{x_{{i - 1}} }}^{{x_{i} }} {a^{j} } (x)\left( {w_{x} (x)} \right)^{j} dx \hfill \\ \,\,\,\,\, + \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {b^{j} (x) - \frac{2}{k}} \right)} w^{j} (x)dx + \int_{{x_{{i - 1}} }}^{{x_{i} }} {\left( {\varpi ^{{j + 1}} (x) + \varpi ^{j} (x)} \right)} dx \hfill \\ \end{gathered}$$

(17)

Using Gauss-Legendre quadrature, we can approximate the integration value in Eq. (17) by setting the integration limits to $-1$ to 1. Now using two-point Gaussian quadrature rule (see [22, 25]), we have

$$\int_{{ - 1}}^{1} {\mathcal{F}} (x)dx \simeq {\mathcal{F}}\left( {\frac{1}{{\sqrt 3 }}} \right) + {\mathcal{F}}\left( { - \frac{1}{{\sqrt 3 }}} \right) + \frac{1}{{135}}{\mathcal{F}}^{{(4)}} (\zeta ),x_{{i - 1}} {\text{ < }}\,\zeta {\text{ < }}\,\,x_{i}$$

(18)

for any function ${\mathcal {F}}$, continuous and differentiable in $[-1,1]$. For function ${\mathcal {F}}(x)$ in an arbitrary interval $\left[ x_{i-1}, x_{i}\right]$, the Gaussian two-point quadrature rule is

$$\int_{{x_{{i - 1}} }}^{{x_{i} }}{\mathcal{F}} (x)dx \simeq \frac{h}{2}\left[ {\mathcal{F}(x_{i} - \kappa ) + {\mathcal{F}}(x_{{i - 1}} + \kappa )} \right],{\text{ }}$$

(19)

where $\kappa =\frac{h}{2}\left( 1-\frac{1}{\sqrt{3}}\right)$. The quadrature rules Eq. (19) are known as Gauss-Legendre rules, and they are based on the roots of Legendre polynomials as their nodes [25]. For many years, the numerical quadrature problem attracted attention because of its many applications, which range from integral equations to the collocation method. An extension of the trapezoidal rule approximation in which the function’s arguments are not predetermined is the two-point Gauss-quadrature rule. The end points $x_{i-1}$ and $x_{i}$ are selected as the quadrature points in the single segment trapezoidal rule for approximating the integral $\int _{x_{i-1}}^{x_{i}}{\mathcal {F}}(x)dx$ whereas in the two point Gaussian-quadrature rule, two points $x_{i-1}$ and $x_{i}$ are selected somewhere between the end points $x_{i-1}$ and $x_{i}$ in order to obtain greater accuracy. Now, using Eq. (19), Eq. (17) can be reduced to:

$$\begin{gathered} - w^{{j + 1}} \left( {x_{i} + \varepsilon _{1} } \right) + w^{{j + 1}} \left( {x_{{i - 1}} + \varepsilon _{1} } \right) + w^{{j + 1}} (x_{i} ) - w^{{j + 1}} (x_{{i - 1}} ) - \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\varepsilon _{2} \left[ {a^{{j + 1}} (x_{i} )w^{{j + 1}} (x_{i} ) - a^{{j + 1}} (x_{{i - 1}} )w^{{j + 1}} (x_{{i - 1}} )} \right] + \hfill \\ \frac{{h\varepsilon _{2} }}{2}\left[ {a^{\prime}\left( {x_{i} - \kappa ,t_{{j + 1}} } \right)w^{{j + 1}} (x_{i} - \kappa ) + a^{\prime}\left( {x_{i} + \kappa ,t_{{j + 1}} } \right)w^{{j + 1}} (x_{{i - 1}} + \kappa )} \right] + \hfill \\ \frac{h}{2}\left[ {\left( {b^{{j + 1}} \left( {x_{i} - k} \right) + \frac{2}{k}} \right)w^{{j + 1}} (x_{i} - \kappa ) + \left( {b^{{j + 1}} \left( {x_{{i - 1}} + k} \right) + \frac{2}{k}} \right)w^{{j + 1}} (x_{{i - 1}} + \kappa )} \right] = \hfill \\ \,\,\,\,\,\,\,\,\,\,\,w^{j} \left( {x_{i} + \varepsilon _{1} } \right) - w^{j} \left( {x_{{i - 1}} + \varepsilon _{1} } \right) - w^{j} (x_{i} ) + w^{j} (x_{{i - 1}} ) + \varepsilon _{2} \left[ {a^{j} (x_{i} )w^{j} (x_{i} ) - } \right. \hfill \\ \left. {a^{j} (x_{{i - 1}} )w^{j} (x_{{i - 1}} )} \right] - \frac{{h\varepsilon _{2} }}{2}a^{\prime}\left( {x_{i} - \kappa ,t_{j} } \right)w^{j} (x_{i} - \kappa ) + a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{j} } \right)w^{j} (x_{{i - 1}} + \kappa ) \hfill \\ - \frac{h}{2}\left[ {\left( {b^{j} \left( {x_{i} - k} \right) - \frac{2}{k}} \right)w^{j} (x_{i} - \kappa ) + \left( {b^{j} \left( {x_{{i - 1}} + k} \right) - \frac{2}{k}} \right)w^{j} (x_{{i - 1}} + \kappa )} \right] + \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\frac{h}{2}\left[ {\varpi ^{{j + 1}} (x_{i} - \kappa ) + \varpi ^{{j + 1}} (x_{{i - 1}} + \kappa ) + \varpi ^{j} (x_{i} - \kappa ) + \varpi ^{j} (x_{{i - 1}} + \kappa )} \right] \hfill \\ \end{gathered}$$

(20)

Since delay term is of order $\varepsilon _{1}$, to tackle the term containing delay, we use Taylor’s approximation of $w^{j+1}\left( x_{i}+\varepsilon _{1}\right) , w^{j+1}\left( x_{i-1}+\varepsilon _{1}\right) , w^{j}\left( x_{i}+\varepsilon _{1}\right) ,w^{j}\left( x_{i-1}+\varepsilon _{1}\right)$,

$w^{j+1}\left( x_{i}-\kappa \right) ,w^{j+1}\left( x_{i-1}+\kappa \right)$, $w^{j}\left( x_{i-1}+\kappa \right)$ and approximating the first derivatives using second order interpolation and using notation $w^{j+1}(x_{i})=W^{j+1}_{i}$, Eq. (20) reduced to:

$$\left\{ \begin{gathered} \left[ { - \frac{{\varepsilon _{1} }}{h} + \varepsilon _{2} a_{{i - 1}}^{{j + 1}} + \frac{{\varepsilon _{2} \kappa }}{4}a^{\prime}\left( {x_{i} - \kappa ,t_{{j + 1}} } \right) + \frac{{\varepsilon _{2} }}{4}\left( {2h - 3\kappa } \right)a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{{j + 1}} } \right)} \right. \hfill \\ \left. { + \frac{\kappa }{4}\left( {b^{{j + 1}} \left( {x_{i} - \kappa } \right) + \frac{2}{k}} \right) + \frac{1}{4}\left( {2h - 3\kappa } \right)\left( {b^{{j + 1}} \left( {x_{{i - 1}} + \kappa } \right) + \frac{2}{k}} \right)} \right]W_{{i - 1}}^{{j + 1}} + \hfill \\ \left[ {\frac{{2\varepsilon _{1} }}{h} - \varepsilon _{2} a_{i}^{{j + 1}} + \frac{{h\varepsilon _{2} }}{2}a^{\prime}\left( {x_{i} - \kappa ,t_{{j + 1}} } \right) + \varepsilon _{2} \kappa a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{{j + 1}} } \right) + } \right. \hfill \\ \left. {\frac{h}{2}\left( {b^{{j + 1}} \left( {x_{i} - \kappa } \right) + \frac{2}{k}} \right) + \kappa \left( {b^{{j + 1}} \left( {x_{{i - 1}} + \kappa } \right) + \frac{2}{k}} \right)} \right]W_{i}^{{j + 1}} + \hfill \\ \left[ { - \frac{{\varepsilon _{1} }}{h} - \frac{{\varepsilon _{2} \kappa }}{4}a^{\prime}\left( {x_{i} - \kappa ,t_{{j + 1}} } \right) - \frac{{\varepsilon _{2} \kappa }}{4}a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{{j + 1}} } \right) - } \right. \hfill \\ \left. {\frac{\kappa }{4}\left( {b^{{j + 1}} \left( {x_{i} - \kappa } \right) + \frac{2}{k}} \right) - \frac{\kappa }{4}\left( {b^{{j + 1}} \left( {x_{{i - 1}} + \kappa } \right) + \frac{2}{k}} \right)} \right]W_{{i + 1}}^{{j + 1}} = \hfill \\ \left[ {\frac{{\varepsilon _{1} }}{h} - \varepsilon _{2} a_{{i - 1}}^{j} - \frac{{\varepsilon _{2} \kappa }}{4}a^{\prime}\left( {x_{i} - \kappa ,t_{j} } \right) - \frac{{\varepsilon _{2} }}{4}\left( {2h - 3\kappa } \right)a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{j} } \right) - } \right. \hfill \\ \left. {\frac{\kappa }{4}\left( {b^{j} \left( {x_{i} - \kappa } \right) - \frac{2}{k}} \right) - \frac{1}{4}\left( {2h - 3\kappa } \right)\left( {b^{j} \left( {x_{{i - 1}} + \kappa } \right) - \frac{2}{k}} \right)} \right]W_{{i - 1}}^{j} + \hfill \\ \left[ { - \frac{{2\varepsilon _{1} }}{h} + \varepsilon _{2} a_{i}^{j} - \frac{{h\varepsilon _{2} }}{2}a^{\prime}\left( {x_{i} - \kappa ,t_{j} } \right) - \varepsilon _{2} \kappa a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{j} } \right) - } \right. \hfill \\ \left. {\frac{h}{2}\left( {b^{j} \left( {x_{i} - \kappa } \right) - \frac{2}{k}} \right) - \kappa \left( {b^{j} \left( {x_{{i - 1}} + \kappa } \right) - \frac{2}{k}} \right)} \right]W_{i}^{j} + \hfill \\ \left[ {\frac{{\varepsilon _{1} }}{h} + \frac{{\varepsilon _{2} \kappa }}{4}a^{\prime}\left( {x_{i} - \kappa ,t_{j} } \right) + \frac{{\varepsilon _{2} \kappa }}{4}a^{\prime}\left( {x_{{i - 1}} + \kappa ,t_{j} } \right) + } \right. \hfill \\ \left. {\frac{\kappa }{4}\left( {b^{j} \left( {x_{i} - \kappa } \right) - \frac{2}{k}} \right) + \frac{\kappa }{4}\left( {b^{j} \left( {x_{{i - 1}} + \kappa } \right) - \frac{2}{k}} \right)} \right]W_{{i + 1}}^{j} + \hfill \\ \frac{h}{2}\left[ {\varpi ^{{j + 1}} (x_{i} - \kappa ) + \varpi ^{{j + 1}} (x_{{i - 1}} + \kappa ) + \varpi ^{j} (x_{i} - \kappa ) + \varpi ^{j} (x_{{i - 1}} + \kappa )} \right]. \hfill \\ \end{gathered} \right.$$

(21)

The stability and convergence of numerical methods can be greatly influenced by the initial condition in the context of TP-SPPs. More robustly, methods that are intended to be uniformly stable-like adaptive meshing or those fitted operators can handle a larger range of initial conditions. In contrast to standard methods, methods such as exponentially fitted schemes or fitted operator methods can achieve uniform convergence, which is less sensitive to the smoothness of the initial condition. To treat the effect of perturbation parameter, exponential fitting factor $(\sigma (\rho ))$ is multiplied Eq. (21) on the term containing $\varepsilon _{1}$ as:

$$\begin{aligned} {\left\{ \begin{array}{ll} \left[ -\dfrac{\varepsilon _{1}\sigma (\rho )}{h}+\varepsilon _{2}a^{j+1}_{i-1}+\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) +\dfrac{\varepsilon _{2}}{4}\left( 2h-3\kappa \right) a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) \right. \\ \left. \hspace{0.5cm}+\frac{\kappa }{4}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) +\frac{1}{4}\left( 2h-3\kappa \right) \left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) \right] W^{j+1}_{i-1}+\\ \left[ \dfrac{2\varepsilon _{1}\sigma (\rho )}{h}-\varepsilon _{2}a^{j+1}_{i}+\dfrac{h\varepsilon _{2}}{2}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) +\varepsilon _{2}\kappa a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) +\right. \\ \left. \hspace{1cm}\frac{h}{2}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) + \kappa \left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) \right] W^{j+1}_{i}+\\ \left[ -\dfrac{\varepsilon _{1}\sigma (\rho )}{h}-\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) -\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) \right. \\ \left. \hspace{0.4cm}-\frac{\kappa }{4}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) -\frac{\kappa }{4}\left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) \right] W^{j+1}_{i+1}=\\ \left[ \dfrac{\varepsilon _{1}\sigma (\rho )}{h}-\varepsilon _{2}a^{j}_{i-1}-\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j}\right) -\dfrac{\varepsilon _{2}}{4}\left( 2h-3\kappa \right) a^{'}\left( x_{i-1}+\kappa ,t_{j}\right) \right. \\ \left. \hspace{0.5cm}-\frac{\kappa }{4}\left( b^{j}\left( x_{i}-\kappa \right) -\dfrac{2}{k}\right) -\frac{1}{4}\left( 2h-3\kappa \right) \left( b^{j}\left( x_{i-1}+\kappa \right) -\dfrac{2}{k}\right) \right] W^{j}_{i-1}+\\ \left[ -\dfrac{2\varepsilon _{1}\sigma (\rho )}{h}+\varepsilon _{2}a^{j}_{i}-\dfrac{h\varepsilon _{2}}{2}a^{'}\left( x_{i}-\kappa ,t_{j}\right) -\varepsilon _{2}\kappa a^{'}\left( x_{i-1}+\kappa ,t_{j}\right) -\right. \\ \left. \hspace{2cm}\frac{h}{2}\left( b^{j}\left( x_{i}-\kappa \right) -\dfrac{2}{k}\right) -\kappa \left( b^{j}\left( x_{i-1}+\kappa \right) -\dfrac{2}{k}\right) \right] W^{j}_{i}+\\ \left[ \dfrac{\varepsilon _{1}\sigma (\rho )}{h}+\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j}\right) +\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i-1}+\kappa ,t_{j}\right) +\right. \\ \left. \hspace{1cm}\frac{\kappa }{4}\left( b^{j}\left( x_{i}-\kappa \right) -\dfrac{2}{k}\right) + \frac{\kappa }{4}\left( b^{j}\left( x_{i-1}+\kappa \right) -\dfrac{2}{k}\right) \right] W^{j}_{i+1}+\\ \hspace{0.4cm}\dfrac{h}{2}\left[ \varpi ^{j+1}(x_{i}-\kappa )+\varpi ^{j+1}(x_{i-1}+\kappa ) +\varpi ^{j}(x_{i}-\kappa )+\varpi ^{j}(x_{i-1}+\kappa )\right] . \end{array}\right. } \end{aligned}$$

(22)

The value of $\sigma (\rho )$ is found in [9] and is given as:

$$\sigma (\rho ) = \left\{ {\begin{array}{*{20}c} {\frac{{\rho \varepsilon _{2} a^{{j + 1}} (0)}}{2}\coth \left( {\frac{{\varepsilon _{2} a^{{j + 1}} (0)\rho }}{2}} \right),} & {{\text{ for}}\;x = 0,} \\ {\frac{{\rho \varepsilon _{2} a^{{j + 1}} (1)}}{2}\coth \left( {\frac{{\varepsilon _{2} a^{{j + 1}} (1)\rho }}{2}} \right),} & {{\text{ for}}\;x = 1,} \\ \end{array} } \right.{\text{ }}$$

(23)

where $\rho =\dfrac{h}{\varepsilon _{1}}$. Eq. (22) can be reformulated as a three-term recurrence relation as:

$$\left\{ {\begin{array}{*{20}l} {\mathbb{L}_{{\varepsilon _{1} ,\varepsilon _{2} }}^{{N,M}} W_{i}^{{j + 1}} = A_{i}^{ - } W_{{i - 1}}^{{j + 1}} + A_{i}^{c} W_{i}^{{j + 1}} + A_{i}^{ + } W_{{i + 1}}^{{j + 1}} = } \hfill \\ {B_{i}^{ - } W_{{i - 1}}^{j} + B_{i}^{c} W_{i}^{j} + B_{i}^{ + } W_{{i + 1}}^{j} + F_{j}^{j} ,} \hfill \\ {W_{i}^{0} = \zeta _{b} (x_{i} ),i = 1,...,N - 1,} \hfill \\ {W_{0}^{{j + 1}} = \zeta _{l} (t_{{j + 1}} ),\quad W_{{N + 1}}^{{j + 1}} = \zeta _{r} (t_{{j + 1}} ),0 \le j \le M - 1.} \hfill \\ \end{array} } \right.{\text{ }}$$

(24)

for $i=0,1,...,N$ and $j=0,1,...,M$ where

$$\begin{aligned} {\left\{ \begin{array}{ll} A^{-}_{i}=-\dfrac{\sigma (\rho )}{\rho }+\varepsilon _{2}a^{j+1}_{i-1}+\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) +\dfrac{\varepsilon _{2}}{4}\left( 2h-3\kappa \right) a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) \\ \hspace{1cm}+\frac{\kappa }{4}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) +\frac{1}{4}\left( 2h-3\kappa \right) \left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) ,\\ A^{c}_{i}=\dfrac{2\sigma (\rho )}{\rho }-\varepsilon _{2}a^{j+1}_{i}+\dfrac{h\varepsilon _{2}}{2}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) +\varepsilon _{2}\kappa a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) +\\ \hspace{2cm}\frac{h}{2}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) +\kappa \left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) ,\\ A^{+}_{i}=-\dfrac{\sigma (\rho )}{\rho }-\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) -\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) -\\ \hspace{2.2cm}\frac{\kappa }{4}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) -\frac{\kappa }{4}\left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) ,\\ B^{-}_{i}=\dfrac{\sigma (\rho )}{\rho }-\varepsilon _{2}a^{j}_{i-1}-\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j}\right) -\dfrac{\varepsilon _{2}}{4}\left( 2h-3\kappa \right) a^{'}\left( x_{i-1}+\kappa ,t_{j}\right) \\ \hspace{1cm}-\frac{\kappa }{4}\left( b^{j}\left( x_{i}-\kappa \right) -\dfrac{2}{k}\right) - \frac{1}{4}\left( 2h-3\kappa \right) \left( b^{j}\left( x_{i-1}+\kappa \right) -\dfrac{2}{k}\right) ,\\ B^{c}_{i}=-\dfrac{2\sigma (\rho )}{\rho }+\varepsilon _{2}a^{j}_{i}-\dfrac{h\varepsilon _{2}}{2}a^{'}\left( x_{i}-\kappa ,t_{j}\right) -\varepsilon _{2}\kappa a^{'}\left( x_{i-1}+\kappa ,t_{j}\right) -\\ \hspace{2.4cm}\frac{h}{2}\left( b^{j}\left( x_{i}-\kappa \right) -\dfrac{2}{k}\right) - a\kappa \left( b^{j}\left( x_{i-1}+\kappa \right) -\dfrac{2}{k}\right) ,\\ B^{+}_{i}=-\dfrac{\sigma (\rho )}{\rho }+\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j}\right) +\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i-1}+\kappa ,t_{j}\right) +\\ \hspace{1.4cm}\frac{\kappa }{4}\left( b^{j}\left( x_{i}-\kappa \right) -\dfrac{2}{k}\right) +\frac{\kappa }{4}\left( b^{j}\left( x_{i-1}+\kappa \right) -\dfrac{2}{k}\right) ,\\ F^{j}_{i}=\dfrac{h}{2}\left[ \varpi ^{j+1}(x_{i}-\kappa )+\varpi ^{j+1}(x_{i-1}+\kappa )+\varpi ^{j}(x_{i}-\kappa )+\varpi ^{j}(x_{i-1}+\kappa )\right] . \end{array}\right. } \end{aligned}$$

(25)

Fitted operator FDM obtained through the two-point Gaussian quadrature rule is the required scheme developed in Eq. (24) to solve the problem in Eq. (1). By ensuring that the initial conditions $\zeta _{b}(x)$ are incorporated into the scheme, these methods offer improved stability and convergence. The matrix inversion method is used to solve Eq. (24). For small mesh sizes, it can be observed from the tridiagonal system of Eqs. (24) and Eqs. (25) that

$$\left| {A_{i}^{c} } \right| \ge \left| {A_{i}^{ - } } \right| + \left| {A_{i}^{ + } } \right|.{\text{ }}$$

(26)

This shows the coefficient matrix associated to ${\mathbb {L}}^{N,M}_{\varepsilon _{1},\varepsilon _{2}}$ is irreducibly diagonally dominant and non-singular. As a result, the matrix is an $M-$matrix, and with the specified boundary conditions, matrix inverse can be used to solve Eq. (24).

Parameter-uniform convergence analysis

We have shown that the continuous solution and its derivatives are bounded. We also evaluated and handled errors resulting from the time variable’s discrete approximation. To guarantee the stability and consistency of developed scheme, we have examined the error estimate in the spatial variable, along with the total discrete scheme.

Lemma 8

(Discrete Maximum Principle). Assume $\Xi ^{j+1}(x)$ be a mesh function satisfying, $\Xi (0,t_{j+1})\ge 0,\Xi (1,t_{j+1})\ge 0$, and ${\mathbb {L}}^{N,M}_{\varepsilon _{1},\varepsilon _{2}}\Xi (x_{i},t_{j+1})\ge 0$ on $\Lambda _{x}^{N}\times \Lambda _{t}^{M}$. Then $\Xi (x_{i},t_{j})\ge 0$ at each point ${\overline{\Lambda }}_{x}^{N}\times {\overline{\Lambda }}_{t}^{M}$, where $\Lambda _{x}^{N}$ and $\Lambda _{t}^{M}$ are discretized domain.

Proof

Let $(x_{\vartheta },t_{j+1})\in \left\{ (x_{i},t_{j+1}):i=1(1)N \right\}$ such that

$$\begin{aligned} \Xi ^{j+1}(x_{\vartheta })=\min _{1 \le i \le N}\Xi ^{j+1}(x_{i}), \end{aligned}$$

and let $\Xi ^{j+1}(x_{\vartheta })<0$. From this $\left( x_{\vartheta },t_{j+1} \right) \notin \left\{ (0,t_{j+1}),(1,t_{j+1})\right\}$, which implies that $(x_{\vartheta },t_{j+1})\in \left\{ (x_{i},t_{j+1}):i=1(1)N-1 \right\}$. Then, from Eq. (24), we have

$$\begin{aligned} {\mathbb {L}}^{N,M}_{\varepsilon _{1},\varepsilon _{2}}\Xi _{\vartheta }^{j+1}=A_{i}^{-}\Xi _{\vartheta -1}^{j+1}+A_{i}^{c}\Xi _{\vartheta }^{j+1}+A_{i}^{+}\Xi _{\vartheta +1}^{j+1} \le 0, \end{aligned}$$

since $\left| A^{-}_{i}\right|>0,\left| A^{c}_{i}\right| >0$ and $\left| A^{+}_{i}\right| >0$, which is a contradiction. As a result, $\Xi ^{j+1}(x_{\vartheta })\ge 0$. Therefore $\Xi ^{j+1}(x_{i})\ge 0, \forall (x_{i},t_{j+1})\in \left\{ (x_{i},t_{j+1}):i=1(1)N-1\right\}$. $\square$

Lemma (8) implies that the discrete operator ${\mathbb {L}}^{N.M}_{\varepsilon _{1},\varepsilon _{2}}$ fulfill to the discrete maximum principle, leading to a monotone coefficient matrix. Moreover, being irreducibly diagonally dominant, it qualifies as an $M-$ matrix. This guarantees the presence of a distinct discrete solution.

Lemma 9

(Estimating Discrete Uniform Stability:) The solution $W_{i}^{j+1}$ of the discrete scheme in Eq. (24) on ${\mathbb {L}}^{N,M}_{\varepsilon _{1}, \varepsilon _{2}}$ gratifies the following estimate:

$$\left\| {W_{i}^{{j + 1}} } \right\| \le \xi ^{{ - 1}} \left\| {\mathbb{L}_{{\varepsilon _{1} ,\varepsilon _{2} }}^{{N,M}} W_{i}^{{j + 1}} } \right\| + \max \left\{ {\left| {\zeta _{l} (t_{{j + 1}} )} \right|,\left| {\zeta _{r} (t_{{j + 1}} )} \right|,\left| {\zeta _{b} (x_{i} )} \right|} \right\},{\text{ }}$$

(27)

for $i=1(1)N-1$ and $0<\xi \le b(x,t)$.

Proof

Let $\varGamma =\xi ^{-1}\left\| {\mathbb {L}}^{N,M}_{\varepsilon _{1},\varepsilon _{2}}W_{i}^{j+1} \right\| +\max \left\{ \left| \zeta _{l}(t_{j+1}) \right| ,\left| \zeta _{r}(t_{j+1}) \right| ,\left| \zeta _{b}(x_{i}) \right| \right\}$ and describe the barrier functions as:

$$\left( {{\triangledown }^{ \pm } } \right)_{i}^{{j + 1}} = \Gamma \pm W_{i}^{{j + 1}} .{\text{ }}$$

(28)

The discrete function $\bigtriangledown ^{\pm }(x_{i},t_{j+1})$ on the boundary and initial conditions are

$$\begin{aligned}&\bigtriangledown ^{\pm }(0,t_{j+1})=\varGamma +\max \left\{ \left| \Upsilon _{l}(t_{j+1}) \right| ,\left| \zeta _{r}(t_{j+1}) \right| ,\left| \zeta _{b}(0) \right| \right\} \pm \zeta _{l}(t_{j+1})\ge 0,\\&\bigtriangledown ^{\pm }(1,t_{j+1})=\varGamma +\max \left\{ \left| \zeta _{l}(t_{j+1}) \right| ,\left| \zeta _{r}(t_{j+1}) \right| ,\left| \zeta _{b}(1) \right| \right\} \pm \Upsilon _{r}(t_{j+1})\ge 0,\\&\bigtriangledown ^{\pm }(x_{i},0)=\varGamma +\max \left\{ \left| \zeta _{l}(0) \right| ,\left| \zeta _{r}(0) \right| ,\left| \zeta _{b}(x_{i},0) \right| \right\} \pm \zeta _{b}(x_{i},0)\ge 0. \end{aligned}$$

Additionally, on the discretized domain, it is given as

$$\begin{aligned} \begin{aligned} {\mathbb {L}}^{N,M}_{\varepsilon _{1},\varepsilon _{2}}\left( \bigtriangledown ^{\pm } \right) ^{j+1}_{i}=\varGamma \pm W^{j+1}_{i}&=A_{i}^{-}\left( \varGamma \pm W_{i-1}^{j+1}\right) +A_{i}^{c}\left( \varGamma \pm W_{i}^{j+1}\right) +A_{i}^{+}\left( \varGamma \pm W_{i+1}^{j+1}\right) ,\\&=\left( A_{i}^{-}+A_{i}^{c}+A_{i}^{+} \right) \varGamma \pm H^{j}_{i}\ge 0, \end{aligned} \end{aligned}$$

where $H^{j}_{i}$ is right hand side of Eq. (24). The proof is right away when Lemma (8) is applied and initial conditions $\zeta _{b}(x)$ are integrated into the scheme. $\square$

The Error of Gaussian quadrature

The definite integral is given by $I_{{\mathcal {F}}}=\int _{x_{i-1}}^{x_{i}}{\mathcal {F}}(x)dx$ and Gaussian quadrature with change of intervals is: $G({\mathcal {F}})=\dfrac{h}{2}\sum _{i=0}^{n}\left[ A_{i}{\mathcal {F}}\left( \dfrac{x_{i}-x_{i-1}}{2}o_{i}+\frac{x_{i}+x_{i-1}}{2} \right) \right]$, where $o_{i}$ is the Gauss nodes on $[-1,1]$. The error of Gaussian quadrature rule is given by [27]:

$$\begin{aligned} E({\mathcal {F}})=I({\mathcal {F}})-G({\mathcal {F}})=\dfrac{h^{6}}{(7!)400}W_{xxxx}(\varphi ),\quad x_{i-1}<\varphi <x_{i}. \end{aligned}$$

To establish a parametric uniform convergence of the discrete scheme of Eq. (24), let $TE^{j+1}(h)$ represent a LTE of the proposed discrete scheme. The LTE is given as:

$$TE^{{j + 1}} (h) = \frac{{h^{2} }}{2}\left[ {\frac{{\varepsilon _{2} \kappa }}{2}a^{\prime}(x_{i} - \kappa ,t_{{j + 1}} ) + \varepsilon _{2} \left( {\frac{{h - 2\kappa }}{2}} \right)a^{\prime}(x_{{i - 1}} + \kappa ,t_{{j + 1}} )} \right] + O(h^{3} ).{\text{ }}$$

(29)

We can rewrite Eq. (24) as follows:

$$\left( { - \frac{{\sigma (\rho )}}{\rho } + \omega _{i} } \right)W_{{i - 1}}^{{j + 1}} + \left( {2\frac{{\sigma (\rho )}}{\rho } + \mu _{i} } \right)W_{i}^{{j + 1}} + \left( { - \frac{{\sigma (\rho )}}{\rho } + z_{i} } \right)W_{{i + 1}}^{{j + 1}} + {\mathcal{K}}_{i} + TE_{i}^{{j + 1}} ,{\text{ }}$$

(30)

where ${\mathcal {K}}_{i}$ is the right hand side of Eq. (24) and

$$\begin{aligned} {\left\{ \begin{array}{ll} \omega _{i}=-\varepsilon _{2}a^{j+1}_{i-1}+\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) +\dfrac{\varepsilon _{2}}{4}\left( 2h-3\kappa \right) a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) +\\ \hspace{1.2cm}\frac{\kappa }{4}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) +\frac{1}{4}\left( 2h-3\kappa \right) \left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) ,\\ \mu _{i}=-\varepsilon _{2}a^{j+1}_{i}+\dfrac{h\varepsilon _{2}}{2}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) +\varepsilon _{2}\kappa a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) +\\ \hspace{1cm}\frac{h}{2}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) +\kappa \left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) ,\\ z_{i}=-\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i}-\kappa ,t_{j+1}\right) -\dfrac{\varepsilon _{2}\kappa }{4}a^{'}\left( x_{i-1}+\kappa ,t_{j+1}\right) -\\ \hspace{1cm}\frac{\kappa }{4}\left( b^{j+1}\left( x_{i}-\kappa \right) +\dfrac{2}{k}\right) - \frac{\kappa }{4}\left( b^{j+1}\left( x_{i-1}+\kappa \right) +\dfrac{2}{k}\right) , \end{array}\right. } \end{aligned}$$

(31)

By imposing the boundary conditions $W^{j+1}(0)=\zeta _{l}(0,t^{j})$ and $W^{j+1}(1)=\zeta _{r}(1,t^{j})$, we obtain system of equations for each time level $t_{j+1}$, in matrix form as:

$$\begin{aligned} \left[ {\mathcal {A}}+{\mathcal {B}}\right] W +{\mathcal {X}}+TE^{j+1}(h)=0, \end{aligned}$$

(32)

where

$$\begin{aligned}&{\mathcal {A}}=trid\left( -\dfrac{\sigma (\rho )}{\rho },2\dfrac{\sigma (\rho )}{\rho },-\dfrac{\sigma (\rho )}{\rho } \right) ,\quad {\mathcal {B}}=trid\left( \omega _{i-1},\mu _{i},z_{i+1} \right) ,\\&{\mathcal {X}}=\left[ {\mathcal {K}}_{1}+\left( \dfrac{\sigma (\rho )}{\rho }+\omega _{1} \right) \zeta _{l}^{j+1}(0),{\mathcal {K}}_{2},...,{\mathcal {K}}_{N-1}+\left( \dfrac{\sigma (\rho )}{\rho }+{\mathcal {K}}_{N-1} \right) \zeta _{r}^{j+1}(1) \right] ,\\&W=\left[ W_{1}^{j+1}, W_{2}^{j+1},...,W_{N-1}^{j+1} \right] ^{T},\quad TE^{j+1}_{i}(h)=\left[ TE^{j+1}_{1},...,TE^{j+1}_{N-1} \right] ^{T}. \end{aligned}$$

Let ${\overline{W}}=\left[ {\overline{W}}_{1}^{j+1}, {\overline{W}}_{2}^{j+1},...,{\overline{W}}_{N-1}^{j+1}\right] ^{T} \cong W$ satisfy the following:

$$\begin{aligned} \left[ {\mathcal {A}}+{\mathcal {B}}\right] {\overline{U}} +{\mathcal {X}}=0, \end{aligned}$$

(33)

Subtracting Eq. (32) from Eq. (33), we have

$$\begin{aligned} \left[ {\mathcal {A}}+{\mathcal {B}}\right] \Psi =TE^{j+1}(h), \end{aligned}$$

(34)

where $\Psi ={\overline{W}}-W={\overline{W}}^{j+1}_{i}-W^{j+1}_{i}=\left[ e_{1},e_{2},...,e_{i-1} \right] ^{T}$ is the discretization error. Since $\varepsilon _{1}$ is small, for a given small h, the matrix $\left[ {\mathcal {A}}+{\mathcal {B}}\right]$ is irreducible and monotone [22]. Hence, there exists $\left[ {\mathcal {A}}+{\mathcal {B}}\right] ^{-1}$ and $\left[ {\mathcal {A}}+{\mathcal {B}}\right] ^{-1}\ge 0$. By manipulating matrix properties and the principles of matrix multiplication and norms, including their inverses, we establish that $\left\| \Psi \right\| \le O(h^{2})$. This implies that the spatial semi-discretization of the proposed method exhibits quadratic convergence.

Theorem 1

Let the solutions to the continuous problem Eq. (1) and the discrete problem Eq. (24) be $w(x_{i},t_{j+1})$ and $W^{j+1}_{i}$ respectively. Next, the fully discrete scheme’s error estimate is provided by

$$\begin{aligned} \max _{i,j}\left| w(x_{i},t_{j+1})-W_{i}^{j+1} \right| \le C\left( k^{2} +h^{2}\right) . \end{aligned}$$

(35)

The proposed method is therefore convergent independent of the perturbation parameters and its rate of convergence is two both in time and space variable.

Proof

The proof is from triangular inequality, Lemma (5) and Lemma (6). $\square$

Numerical examples and results

In this section, the proposed method is implemented on two test problems. To support the theoretical analysis, we have provided the tabular results for errors and order of convergence. Since the exact solution is unknown, the maximum point-wise error $E_{\varepsilon _{1},\varepsilon _{2}}^{N,M}$ and rate of convergence ${\overline{R}}_{\varepsilon _{1},\epsilon _{2}}^{N,M}$ of the proposed scheme are computed using the double mesh principle with given initial and boundary conditions. The maximum point-wise absolute error and the convergence rate are defined as follows:

$$\begin{aligned} E_{\varepsilon _{1},\varepsilon _{2}}^{N,M}=\max _{0\le i \le N,0\le j \le M}\left| W_{i,j}^{N,M}-W_{i,j}^{2N,2M} \right| , {\overline{R}}_{\varepsilon _{1},\varepsilon _{2}}^{N,M}=\dfrac{\log \left( E_{\varepsilon _{1},\varepsilon _{2}}^{N,M}\right) -\log \left( E_{\varepsilon _{1},\varepsilon _{2}}^{2N,2M} \right) }{\log (2)}, \end{aligned}$$

where $W_{i,j}^{N,M}$ and $W_{i,j}^{2N,2M}$ are numerical solutions computed on the mesh $N\times M$ and $2N\times 2M$, respectively. Additionally, the uniform rate of convergence ${\overline{R}}^{N,M}$ and uniform maximum point-wise error $E^{N,M}$ are given by the following formula:

$$\begin{aligned} E^{N,M}=\max _{\varepsilon _{1},\epsilon _{2}}E_{\varepsilon _{1},\varepsilon _{2}}^{N,M},\quad {\overline{R}}^{N,M}=\dfrac{\log \left( E^{N,M}\right) -\log \left( E^{2N,2M} \right) }{\log (2)}. \end{aligned}$$

Example 1

Consider the problems in [9, 13]

$$\begin{aligned} {\left\{ \begin{array}{ll} \dfrac{\partial w}{\partial t}-\varepsilon _{1} \dfrac{\partial ^{2} w}{\partial x^{2}}-\varepsilon _{2}(1+\exp (x))\dfrac{\partial w}{\partial x}+(1+x^{5})w=-10\exp (t^{2})x^{2}(1-x^{2}),\\ w(x,0)=0,w(0,t)=0=w(1,t). \end{array}\right. } \end{aligned}$$

(36)

Example 2

Consider the problems in [15, 28]

$$\begin{aligned} {\left\{ \begin{array}{ll} \dfrac{\partial w}{\partial t}-\varepsilon _{1} \dfrac{\partial ^{2} w}{\partial x^{2}}-\varepsilon _{2}(1+x-x^{2}+t^{2})\dfrac{\partial w}{\partial x}+(1+5xt)w=(x^{2}-x)(e^{t}-1),\\ w(x,0)=0,w(0,t)=0=w(1,t). \end{array}\right. } \end{aligned}$$

(37)

Table 1 Computed $E_{\varepsilon _{1},\varepsilon _{2}}^{N,M},E^{N,M},{\overline{R}}_{\varepsilon _{2},\varepsilon _{2}}^{N,M}$ and ${\overline{R}}^{N,M}$ for Example (1) at $\varepsilon _{2}=10^{-6}$

Table 2 Computed maximum point-wise error and rate of convergence for Example(1) at $\varepsilon _{2}=10^{-6}$

Table 3 Computed maximum point-wise error for Example (1) for different values of $\varepsilon _{1}$ and $\varepsilon _{2}$ by taking $N=64$ and $M=32$

Table 4 Comparison of $E^{N,M}_{\varepsilon _{1},\varepsilon _{2}}, E^{N,M},{\overline{R}}^{N,M}$ for Example (1) with [9, 13]

Table 5 Computed maximum point-wise error and rate of convergence for Example(1) at $\varepsilon _{1}=10^{-6}$ and varying values of $\varepsilon _{1}$

Table 6 Maximum point-wise error & rate of convergence for Example (2) at $\varepsilon _{2}=10^{-7}$

Table 7 Computed maximum point-wise error for Example(2) for different values of $\varepsilon _{1}$ and $\varepsilon _{2}$ by taking $N=128$ and $M=32$

Table 8 Comparison of $E^{N,M}_{\varepsilon _{1},\varepsilon _{2}}, E^{N,M}$ for Example(2) with [15, 28]

Table 9 Computed $E^{N,M}_{\varepsilon _{1},\varepsilon _{2}},E^{N,M},{\overline{R}}^{N,M}$ for Example(2) at $\varepsilon _{1}=10^{-9}$ for $M=N$

Fig. 1

The numerical solution profile: (a) for $N=M=64, \varepsilon _{1}=10^{-2},\varepsilon _{2}=10^{-8}$ and (b) for $N=M=128, \varepsilon _{1}=10^{-12},\varepsilon _{2}=10^{-6}$ for Example(1)

Fig. 2

The numerical solution profile: (a) for $N=M=64, \varepsilon _{1}=10^{-4},\varepsilon _{2}=10^{-2}$ and (b) for $N=M=128, \varepsilon _{1}=10^{-8},\varepsilon _{2}=10^{-6}$ for Example (2)

Fig. 3

Log-log plot for Example (1): (a) for Table 1) and (b) for Table 5

Fig. 4

Log-log plot for Example (2): (a) for Table 7 and (b) for Table 9

The computed maximum point-wise error $E_{\varepsilon _{1},\varepsilon _{2}}^{N,M}$ and rate of convergence ${\overline{R}}_{\varepsilon _{1},\varepsilon _{2}}^{N,M}$ for Examples (1) and  (2) are presented in Tables 1, 2, 3, 4, 5, 6, 7, 8 and 9 with different values of $\varepsilon _{1},\varepsilon _{2}, N$ and M. Tables (4), (6) and (8) compare numerical results to those given in [9, 13, 15, 28], by fixing or varying $\varepsilon _{1}$ or $\varepsilon _{2}$. This comparison indicates that, once compared to studies found in the literature, the suggested approach provides a more accurate solution. Consequently, tabulated numerical results confirm that our theoretical estimates in Lemma (6) and Theorem (1) can be achieved by the proposed scheme, i.e., two-point Gaussian quadrature rule. As the number of mesh points increases, as we observe in Tables 1, 2, 5, 6, 9) errors decrease, and as we approach downward in each column, errors become constant. These results confirm parameter uniformity and second-order accuracy in the relevant cases. All tables reveal that convergence is independent of perturbation parameters, and that the maximum absolute error decreases as the number of mesh points increases. Tabulated results in (Tables 3, 7) indicats that for a fixed value of N, M and $(\varepsilon _{1},\varepsilon _{2})$, maximum point-wise errors going to stabilized as parameter $\varepsilon _{1}$ or $\varepsilon _{2}$ approaches to zero respectively. The numerical solution surface plots for the test problems (1) and (2) are presented in Figs. 1 and 2 respectively with different values of $N,M,\varepsilon _{1}$ and $\varepsilon _{2}$. The problem displays physical behavior of the solution, i.e., left and right boundary layers based on the size of $\varepsilon _{1}$ or $\varepsilon _{2}$, as a result of the small parameters. In addition, Figs. 3 and 4 show the log-log plot of the maximum point-wise errors for Examples (1) and  (2) respectively. Figure 3a depict the log-log plot of maximum point-wise errors for Example (1) taking the values in Table 1. From these plots one can deduce that the numerical scheme converges $(\varepsilon _{1},\varepsilon _{2})-$uniformly as the perturbation paramters goes very small.

Conclusion and future works

A two-parameter singularly perturbed parabolic convection-diffusion problem is considered via the Crank-Nicolson method for the discretization of time derivative, whereas in the discretization of the spatial variable, a two-point Gaussian quadrature formula is used by introducing a fitting factor. In comparison to more conventional numerical methods, the two-point Gaussian quadrature method with an exponential fitting factor offers a new and practical method for resolving problems under consideration. The method’s novelty resides in its independence from asymptotic expansions and fitted meshes. The presence of two small parameters causes the problem to exhibit left and right boundary layers with different width depending on the size of $\varepsilon _{1}$ or $\varepsilon _{2}$. The method is demonstrated to have parameter uniform second order convergence in both time and space. The proposed scheme’s performance is investigated by comparing the results, and it is discovered that the accuracy of the numerical results is comparable to or better than that of existing difference schemes [9, 13, 15, 28], as verified both theoretically and numerically. To further illustrate the versatility and applicability of the suggested approach, we point out that the analysis method utilized in this work can be extended to TP-SPPPs with robin boundary conditions.