An almost second order uniformly convergent method for a two-parameter singularly perturbed problem with a discontinuous convection coefficient and source term

1.   Introduction

Singularly perturbed problems (SPPs) frequently appear in several fields of applied mathematics, such as fluid dynamics ^[22], chemical reactor theory ^[26], etc. These problems are characterized by the presence of a small positive parameter multiplying the highest-order derivative, but other small parameters may appear affecting other terms. Depending on the parameter location and the smoothness of the coefficients of the SPP, the solution shows steep gradients on the boundary and interior parts of the domain, which are characterized as boundary and interior layers. The layer occurrence causes several hurdles (which include huge computational costs in the case of uniform meshes) in the numerical analysis, which is the reason to consider an adaptive mesh. The simplest adaptive mesh in this context is due to Shishkin ^[9], who proposed a piecewise uniform mesh to capture the layer behavior. In contemporary literature, numerical methods for SPPs (involving only a diffusion parameter) with smooth data can be seen in ^[6,9,20], while for non-smooth data, the works ^[1,7,8,10] considered Shishkin meshes.

Two-parameter problems involving convection ($\varepsilon_c$) and diffusion ($\varepsilon_d$) parameters extend the convection and reaction-dominated models. In recent years, several higher-order accurate boundary layer resolving numerical methods based on hybrid schemes, cubic spline schemes, asymptotic expansion methods, etc., have been presented in ^{[11,12,15,24,25,27]} for two-parameter problems with smooth data. The non-smooth data produce interior layers in addition to the boundary layers, whose sharpness depends on the sign of the convection coefficient and has been investigated by various researchers for both singularly perturbed ordinary and partial differential equations (for instance, see ^{[2,3,4,17,19,21]}). The works ^[7,10,14] (involving only a diffusion parameter) are devoted to the numerical analysis of interior layers due to the presence of a discontinuous convection coefficient. A rigorous analysis of the effect of all possible subclasses of discontinuous convection coefficients can be found in Riordan ^[14,18]. These above works motivated us to develop a higher-order numerical approximation for a two-parameter singularly perturbed problem where some of the coefficients have jump discontinuities.

Motivated by the studies in ^[14,18], we consider the following two-parameter singularly perturbed problem with a discontinuous convection coefficient and a discontinuous source term:

Here $\varepsilon_d$ and $\varepsilon_c$ are known as singular perturbation parameters, where $0 < \varepsilon_d \ll 1$, $0\leq \varepsilon_c\leq 1$. For simplicity, we consider the domain as $\overline{\Gamma} = [0, 1]$, with $\Gamma = (0, 1)$, $\Gamma^- = (0, d)$ and $\Gamma^+ = (d, 1).$ Here $b(x)$ is assumed to be a sufficiently smooth function in $\overline{\Gamma}$ satisfying $b(x)\geq \beta > 0$ and $a(x)$, $f(x)$ are sufficiently smooth in $(\Gamma^- \cup \Gamma^+)\cup \{0, 1\}$. Also, $a(x)$ and $f(x)$ have a jump discontinuity at $d \in \Gamma$, where the jump of $\omega(x)$ at $x = d$ is denoted as $[\omega](d) = \omega(d^+)-\omega(d^-)$. These assumptions ensure that the SPP (1.1)-(1.2) has a solution $u(x)\in C^0(\overline{\Gamma})\cap C^1(\Gamma)\cap C^2(\Gamma^- \cup \Gamma^+)$.

Note that $\varepsilon_c = 1$ reduces the general problem in (1.1) to a convection-diffusion problem ^[7], and $\varepsilon_c = 0$ reduces (1.1) to a reaction-diffusion problem ^[8]. The nature of the solution behaves differently with respect to the ratio of the parameters $\varepsilon_d/\varepsilon_c^2.$ The solution behaves similarly to a dissipative case, when $\varepsilon_d / \varepsilon_c^2 \rightarrow 0$ as $\varepsilon_c \rightarrow 0$, and it acts as in the dispersive case when $\varepsilon_c^2/\varepsilon_d \rightarrow0$ as $\varepsilon_d \rightarrow 0$ ^[16]. Hence, we consider the following two cases for the numerical analysis of (1.1):

Case (i): $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$,

Case (ii): $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$,

where $\gamma = \min_{\overline{\Gamma}} \left\{ \dfrac{b(x)}{\alpha(x)}\right\}$, where $\alpha(x) = \alpha_1, \; x < d$, and $\alpha(x) = \alpha_2, \; x > d$.

Apart from the assumptions considered in (1.3), we have also noted the case when the sign of the convection coefficient is reversed in $\Gamma^-$ and $\Gamma^+$ (i.e., $a(x)\geq \alpha_1 > 0, \; x \in \Gamma^-$ and $a(x) \leq -\alpha_2 < 0, \; x \in \Gamma^+$) as a remark.

For two-parameter singularly perturbed problems with smooth data, Riordan et al. ^[15] analyzed the upwind scheme on a Shishkin mesh and obtained a uniform accuracy of $O(N^{-1} \ln^2 N)$, where $N$ defines the number of partitions in the domain. Later, Gracia et al. ^[11] established a higher-order numerical method for this problem, which is based on the combination of upwind, central difference, and mid-point schemes in various partitions of the domain. Their numerical method provided parameter-uniform convergence of accuracy $O(N^{-2} \ln^3 N)$, if $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, and $O(N^{-2} \ln^2 N)$, if $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$ on a Shishkin mesh. In ^[23], the authors introduced a discontinuity in the source term for this problem and obtained $O(N^{-1} \ln^2 N)$, for $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ and $O(N^{-1} \ln^3 N)$, for $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$ on a Shishkin mesh. Several other adaptive meshes, based on arc-length equidistribution ^[13] and curvature-based monitor functions ^[5], have also been considered to improve uniform accuracy up to the first-order.

The above articles motivated us to develop a higher-order numerical analysis for two-parameter problems with a discontinuous convection coefficient and source term. The discontinuous data make the numerical analysis different as it gives rise to the interior layer in addition to the boundary layer. We consider the error analysis on the Shishkin mesh and show that the error is independent of the convection and diffusion parameters.

Throughout this article, we denote $C$ as a generic positive constant independent of the number of nodal points and the perturbation parameters $\varepsilon_d$, $\varepsilon_c$. The convergence is estimated in the infinity norm, which is denoted as $\|u\|_{\Omega} = \max_{x\in \Omega}|u(x)|$ for a function $u(x)$ defined on a general domain $\Omega$. We also write $\|.\| = \|.\|_{\Omega}$, if the norm and domain are obvious. Accordingly, the corresponding discrete norm is denoted as $\|.\| = \|.\|_{{\Omega}^N}$. Without loss of generality, we assume that the number of mesh intervals $N$ is divisible by $2$.

The paper is arranged as follows: In Section 2, we note the existence of the solution and derive a minimum principle for (1.1)-(1.2), from which it follows the stability of the solution $u(x)$. Some estimates of the solution and its derivatives are also stated here. Section 3 presents a discrete problem based on a hybrid finite difference scheme corresponding to the continuous problem. A decomposition of the discrete solution is introduced in Section 4, which helps us evaluate the truncation error estimate. In Section 5, we have shown that this estimate provides a higher-order $\varepsilon_d$-$\varepsilon_c$ uniform numerical approximation in the discrete maximum norm. Numerical examples, given in Section 6, validate the theoretical findings. In the end, we draw a conclusion by highlighting the major contribution of the paper.

2.   Derivative bounds of the continuous solution

In this section, we consider a few analytical properties of the solution of (1.1)-(1.2). The solution is decomposed into regular and singular components, which describe the solution behavior at the outer and inner regions of the boundary and interior layers, respectively. This decomposition will be used in the subsequent sections to obtain a parameter-uniform error estimate. We begin this section with an existence theorem.

Theorem 2.1. The SPP (1.1)-(1.2) has a solution $u(x)$, which belongs to the class $C^0(\overline\Gamma)\cap C^1(\Gamma)\cap C^2(\Gamma^- \cup\Gamma^+)$.

Proof. By using the constructive method presented in ^[8,10], a solution of the SPP (1.1)-(1.2) can be obtained. □

The operator $L$ at (1.1) satisfies the following minimum principle on $\overline{\Gamma}$.

Lemma 2.1. If a function $u(x) \in C^0(\overline{\Gamma}) \cap C^2(\Gamma^-\cup \Gamma^+)$ satisfies $u(0)\geq 0$, $u(1)\geq 0$, $Lu(x) \leq 0$ for all $x \in (\Gamma^- \cup \Gamma^+)$ and $[u'](d) \leq 0$, then $u(x) \geq 0 \; \forall\; x \in \overline{\Gamma}$.

Proof. For the proof and the existence of the solution of the SPP (1.1)-(1.2), the reader is referred to ^[18]. □

As a consequence, we get the following stability estimate:

Lemma 2.2. Let $u(x)$ be a solution of (1.1)-(1.2), then

One can also obtain the following bounds for the solution derivatives when $a(x)$ and $f(x)$ have a jump discontinuity at $x = d$ (see ^[18]).

Lemma 2.3. Let $u(x)$ be the solution of problem (1.1)-(1.2), where $|u(0)|\leq C$ and $|u(1)| \leq C$. Then, for all $0 \leq k \leq 4$, it holds

Proof. Let us begin the discussion on the domain $\Omega^-$. Consider any point $x \in (0, d)$ and a neighborhood $N_p = (p, p+r)$, where $r$ is positive, such that $x \in N_p \subset (0, d)$. Since $u$ is differentiable in $N_p$, the Mean Value Theorem implies that there exists $q \in N_p$ such that $u'(q) = \frac{u(p+r)-u(p)}{r}$. Now, we obtain

We can define $u'(x)$ as

Further, the result for the first derivative can be deduced assuming that $x-q \leq r$ and $r = \sqrt{\varepsilon_d}$, from which we obtain

The bounds for the second derivative will be

Now, we differentiate (1.1) to obtain the third derivative bounds, which result in

Finally, we derive the bounds for the fourth derivative as follows:

After simplifying it, we can write the required result as

The proof for the domain $\Omega^+$ can be obtained using similar arguments. □

Before going into further details about the decomposition of $u(x)$ into regular and singular components, we will consider the following: Let $F(x)$ be a smooth function in $(\Gamma^-\cup\Gamma^+)$, such that $F(x)$ and its derivatives have a jump discontinuity at $d\in \Gamma$. Consider $u(x) \in {\mathcal C}^1(\Gamma)\cap {\mathcal C}^2(\Gamma^-\cup\Gamma^+)$, such that

It can be proven that problem (2.1) has a unique solution ^[7]. Let

and further let $u_l^*(x)$ be the solution of

Similarly, one can define $u_r^*(x)$ on the interval $[d, 1]$. Now

To establish a sharper bound on the error analysis, the solution $u(x)$ is decomposed into a regular component $v(x)$ and a singular component $w(x)$ such that

and $w(x) = w_l(x)+w_r(x)$ where

We now define the regular and singular components as the solutions to the following problems, respectively.

and

Now, let us consider the case (i): $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$.

We first define $v_0(x), \; v_1(x)$ and $v_2(x)$ as the solutions to the following problems:

Choose $v_3(x) \in C^0(\overline{\Gamma}) \cap C^1(\Gamma) \cap C^2(\Gamma^-\cup \Gamma^+)$, such that

Adopting the procedure from ^[11,18], we obtain the upper bounds of the derivatives of the regular and the singular components given in the following Lemmas 2.4 and 2.5.

Lemma 2.4. When $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, the regular component $v(x)$ and its derivatives satisfy the following bounds:

Lemma 2.5. When $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, the singular components $w_l(x)$ and $w_r(x)$ and their derivatives satisfy the bounds

where

Now consider the case (ii): $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$.

Let $v_0(x), \; v_1(x)$, and $v_2(x)$ be the solutions of the following problems, respectively:

Choose $v_3(x) \in C^0(\overline{\Gamma}) \cap C^1(\Gamma) \cap C^2(\Gamma^-\cup \Gamma^+)$ such that

Using the reasoning given in ^[11,18], the following Lemmas 2.6 and 2.7 can be proved for the case $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$.

Lemma 2.6. When $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, the regular component $v(x)$ satisfies the following bounds

Lemma 2.7. When $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$ the singular components $w_l(x)$ and $w_r(x)$ satisfy the bounds

where

It can be verified that $v(x)+w_l(x)+w_r(x)$ satisfies the problem (1.1)-(1.2). Therefore, the unique solution to the problem is

3.   Discretization of the problem using a Hybrid scheme

In this section, we introduce a difference scheme to discretize the continuous problem (1.1)-(1.2). The discrete problem combines the standard upwind, mid-point, central difference, and five-point difference schemes. The five-point difference scheme is applied at the point of discontinuity to ensure higher-order (in this case, second-order) accuracy. This scheme is reduced to a three-point structure to preserve the monotonicity property with almost second-order accuracy. The numerical scheme, obtained by combining all these operators, maintains the monotonicity property.

The discrete problem will be defined on an a priori adaptive piecewise uniform mesh, which is dense inside the boundary and interior layer regions. To construct this mesh, we first divide the domain $\overline{\Gamma}$ into six subintervals:

for some $\tau_1, \, \tau_2, \tau_3$, and $\tau_4$. The mesh points are denoted as $\overline{\Gamma}^N = \{x_i\}_0^N$, where $x_{N/2}$ denotes the point of discontinuity, $x_{N/2} = d.$ On these mesh points, we define the discrete solution as $U_i.$ The transition parameters $\tau_1, \; \tau_2, \; \tau_3$, and $\tau_4$ in $\overline{\Gamma}$ are chosen as follows:

where $\theta_1$ and $\theta_2$ are defined in the previous section. Now we construct a uniform mesh on each of the subintervals $[0, \tau_1], \; [d-\tau_2, d], \; [d, d+\tau_3]$ and $[1-\tau_4, 1],$ so that each one contains $N/8+1$ uniform mesh points, and the subintervals $[\tau_1, d-\tau_2]$ and $[d+\tau_3, 1-\tau_4]$ contain $N/4+1$ uniform mesh points respectively. The mesh sizes in each of the subintervals from left to right of $\overline{\Gamma}$ are denoted as $h_1 = 8\tau_1/N, \; h_2 = 4(d-\tau_2-\tau_1)/N, \; h_3 = 8\tau_2/N, \; h_4 = 8\tau_3/N, \; h_5 = 4(1-d-\tau_3-\tau_4)/N, \; \mbox{and}\; h_6 = 8\tau_4/N.$ On the above adaptive mesh $\overline{\Gamma}^N,$ we discretize the BVP (1.1)-(1.2) as

The discretization (3.2), based on the upwind scheme, is almost first-order accurate ^[23] (see the numerical Section 6). Hence, our goal is to improve the order of accuracy. Now, we construct an almost second-order accurate hybrid scheme for (1.1)-(1.2) by combining the following central, upwind, and mid-point difference schemes with a five-point scheme:

where, $D^0 U_i = \dfrac{U_{i+1} - U_{i-1}}{h_i + h_{i+1}}, \; D^+ U_i = \dfrac{U_{i+1} - U_i}{h_{i+1}}$, $D^- U_i = \dfrac{U_{i} - U_{i-1}}{h_{i}}$, $\delta^2 U_i = \dfrac{1}{\overline{h}_i} \left(D^+ U_i - D^- U_i\right)$, $\overline{z} _i = \dfrac{z_i + z_{i+1}}{2}$ and $D^* = \begin{cases} D^-, & \mbox{if}\; i < N/2, \\ D^+, & \mbox{if}\; i > N/2. \end{cases}$

At $x_{N/2} = d$, we use a five-point difference scheme by combining the second-order accurate one-sided forward difference approximation $u'(x)\approx (-3 U(x) + 4 U(x+h_3) - U(x+2h_3))/2h_3$ with the backward difference approximation $u'(x)\approx (3 U(x) - 4 U(x-h_4) + U(x-2h_4))/2h_4$. At the point of discontinuity, we define the scheme as:

where $h = \max\{h_3, h_4\}$.

The finite difference scheme, which involves the central, upwind, mid-point, and five-band difference schemes on the piecewise uniform mesh, will be written as

where $L^N$ is defined as follows:

When $\sqrt{\alpha}\varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$,

and for $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$

At the transition points $\tau_1$ and $(d+\tau_3),$ the scheme is defined as

At the transition points $(d-\tau_2)$ and $(1-\tau_4),$ we define the scheme as

Now, we define the discrete scheme as

The matrix associated with (3.4) does not satisfy the M-matrix condition at the point of discontinuity $x_i = d$. But, without loss of generality, we can convert this five-point difference scheme into a three-point difference scheme (say, $L^N_{T}U_i$) by estimating $U_{N/2-2}, \; U_{N/2+2}$ from $L_c^N U_i$ so that the new equations have the monotonicity property. To do this, we take

Now we replace the above expressions of $U _{N/2-2}, \; U _{N/2+2}$ at the five-point difference scheme ($L^N_t U_{N/2}$) to construct a three-point scheme ($L^N_{T} U_{N/2}$) that preserves the monotonicity property and leads to a higher-order accuracy at the point of discontinuity, as given by

So, the reformulated discrete operator (say $L_*^N U_i$) of (3.4) can be written as

where

and

The entries of the stiffness matrix corresponding to $L_*^N$ are given by

Now, we state a few conditions that are used to preserve the monotonic properties of the discrete problem (3.5). In $(0, \tau_1)$ and $(d, d+\tau_3),$ note that

For $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ in $(d-\tau_2, d)$ and $(1-\tau_4, 1),$ we get

For $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$ in $(d-\tau_2, d)$ and $(1-\tau_4, 1),$ we get

At $x_i = d$, we have

We note that in order to guarantee that the operator $L^N_*$ is monotone, it is necessary to impose the following assumption:

which will be evident in the proof of the following lemma. Thus, the discrete problem is

The following lemma shows that the discrete operator $L^N_*$ has stability properties analogous to those of the continuous operator $L$.

Lemma 3.1. If a mesh function $U_i$ satisfies $U_0 \geq 0$, $U_N \geq 0$, $L^N_* U_i \leq 0$ for all $i = 1, 2, \ldots, N-1$ and $L^N_T U_{N/2} \leq 0$, then $U_i \geq 0$ for all $i = 0, 1, \ldots, N$.

Proof. The system $L^N_* U_i = Q^N_* f_i$ for all $1\leq i \leq N-1$ is a linear system of $N-1$ equations. We show that the corresponding matrix is diagonally dominant. To show this, it is sufficient to check that the conditions

are fulfilled for the operators defined in (3.11).

For $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, the central difference operator $(L^N_c)$ is used in all the subintervals of $\overline{\Gamma}^N$. Therefore, $s^-_i > 0$ is guaranteed in the subintervals $(d-\tau_2, d)$ and $(d, d+\tau_3)$ by the definition of $h_3$ and $h_4$. In the remaining regions, the $L^N_c$ operator is used when $\varepsilon_c h_1 \|a\|, \varepsilon_c h_2 \|a\|, \varepsilon_c h_5 \|a\|, \varepsilon_c h_6 \|a\| < 2\varepsilon_d$, and thus, $s^-_i > 0$ is satisfied from (3.7), (3.8), and (3.10).

For $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, the central difference operator is applied in $(d-\tau_2, d)$ and $(d, d+\tau_3).$ Here, $h_3$ and $h_4$ are given by $s^-_i > 0$ and $s^+_i > 0$. The mid-point operator is applied to the layer region $(0, \tau_1)$. Here, the inequality$h_i s^-_i = \dfrac{\varepsilon_d}{\overline{h}_i} - \varepsilon_c a_i > 0$ is guaranteed since $\varepsilon_c h_1 \|a\| < \dfrac{\varepsilon_d}{2}$ and $s^+_i = \dfrac{\varepsilon_d}{h_{i+1}\overline{h}_i} - \dfrac{b_{i+1}}{2} > 0$ is satisfied if $\|b\| h_1^2 < \dfrac{\varepsilon_d}{4}$. The inequalities (3.6) and (3.10) make it obvious. The sign pattern of $s^+_i+s^c_i+s^-_i < 0$ is direct using the condition $s^c_i = -s^+_i-s^-_i-\overline{b}_i$. When the mid-point operator is used for the interval $(1-\tau_4, 1)$ the condition (3.12) is satisfied if $s_i^+ = \dfrac{\varepsilon_c \overline{a}_i}{h_{i+1}} - \dfrac{b_i}{2} > 0.$ This is guaranteed since $\|b\| h_6 < \varepsilon_c \alpha/2$.

In the coarse mesh region $(\tau_1, d-\tau_2)$ and $(d+\tau_3, 1-\tau_4)$, $s^+_i > 0$ also follows from the mid-point operator since $2\|b\|h_2 < \varepsilon_c\alpha$ and $2\|b\| h_5 < \varepsilon_c \alpha$ on each of the intervals. If $2\|b\|h_2, 2\|b\|h_5 \geq \varepsilon_c \alpha$, the upwind operator $L^N_u$ is used to preserve $s^+_i > 0$.

At the point of discontinuity, $x_{N/2} = d$, where $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, $s_{N/2}^+ = 2\varepsilon_d + 2 h_4 \varepsilon_c a_{N/2} - b_{N/2} h_4^2 > 0$ is guaranteed since $\|b\| h_4 \leq 2 \varepsilon_c \alpha$. For $s_{N/2}^- = 2\varepsilon_d - b h_3^2 -2 \varepsilon_c h_3 a_{N/2} > 0$ is assured with $\|b\| h_3^2 < \varepsilon_d /4$ and $\varepsilon_c h_3 \|a\| < \varepsilon_d/2$. Again, $s_{N/2}^c = -s^-_{N/2} - s^+_{N/2} - 2b_{N/2-1} \left[\dfrac{h_3^2}{2\varepsilon_d + \varepsilon_c a_{N/2-1} h_3} + \dfrac{h_4^2}{2\varepsilon_d + \varepsilon_c a_{N/2+1} h_4}\right]$ follows, $s_{N/2}^c < 0$. For $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, we have $s_{N/2}^+ > 0$ and $s_{N/2}^- > 0$ since $\|b\|h_4^2 < \varepsilon_d/4$ and $2\|b\| h_3 < \varepsilon_c \alpha$. Similarly, we can show $s_{N/2}^c < 0$.

Hence, combining all the above relations defined at various mesh points, it can be noted that the matrix corresponding to $L^N_*$ is an M-matrix, and hence it satisfies the discrete minimum principle. □

Lemma 3.2. If $U_i$ is the solution to (3.11), then

where $i = 0, 1, 2, \ldots, N$ and $\theta = \min \left\{\alpha_1/d, \alpha_2/(1-d) \right\}$.

Proof. Let us define the mesh functions

where $M = \max\left\{ |U_0|, |U_N|\right\} + \dfrac{1}{\theta} \|Q^N_* f_i\|$.

It is obvious, $\Theta^{\pm}_0 \geq 0$, $\Theta^{\pm}_N \geq 0$ and also

Similarly, $L_m^N \Theta^{\pm}_i \leq 0$ and $L_u^N \Theta^{\pm}_i \leq 0$ for the above values of $i$.

At the point of discontinuity, it is

Applying Lemma 3.1, it follows that

which leads to the desired bound on $U_i$. □

4.   Truncation error analysis

We address the error analysis in this section. The nodal error is denoted by $e_i = U_i-u(x_i)$. To bound the nodal error $|e_i|,$ we first decompose the discrete solution (in a similar manner, as was done with the continuous solution) of (3.11) as $U_i = V_i+W_{l, i}+W_{r, i}$. The discrete regular component ($V_i$) and singular components ($W_{l, i}$ and $W_{r, i}$) are again decomposed to obtain sharper bounds. Let us define the mesh functions $V^-_i$ and $V^+_i$, which approximate $V_i$ on either side of the point $x_i = d$. In addition, we construct the mesh functions $W_{l, i}^-, \; W_{l, i}^+$ and $W_{r, i}^-, \, W_{r, i}^+$ to approximate $W_{l, i}$ and $W_{r, i}$ on the left and right sides of $x_i = d$. Here, $W_{l, i}^-$ and $W_{r, i}^-$ correspond to the left boundary layer and right interior layer, respectively. Similarly, $W_{l, i}^+$ and $W_{r, i}^+$ are the solutions of the left interior layer and right boundary layer. These mesh functions are useful to show the convergence of the nodal error $|e_i|$ inside and outside the layers.

Let the regular components $V^-_i$, $V^+_i$ are the solutions to the following discrete problems:

Similarly, the discrete problem corresponding to the left singular components $W^-_{l, i}$ and $W^+_{l, i}$ is defined as follows:

The corresponding right singular components $W_{r, i}^-$ and $W_{r, i}^+$ can be described in a similar way.

Consequently, the solution of the discrete problem (3.11) can be written as

The following lemma provides the bounds for the singular components.

Lemma 4.1. The bounds of the left and right singular components $W_{l, i}^-, \; W_{l, i}^+, \; W_{r, i}^-$ and $W_{r, i}^+$ are as follows:

where $k_1 = i$, $k_2 = N/2+i$, and $C$ is a positive constant independent of $\varepsilon_c, \, \varepsilon_d$.

Proof. Define the barrier functions

where

and also

Now we will prove that $L^N_* \psi_{l, i}^- \leq 0, \; L^N_* \psi_{r, i}^- \leq 0$. Applying the discrete operator (3.11) on $\psi_{l, i}^-,$ we have

Again, for the central difference operator, it follows that

Applying (4.1) in the above equation for both the cases $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ and $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, we obtain

In a similar way, the upwind and mid-point difference schemes on $\psi_{l, i}^-$ lead to

and

Combining all the above inequalities, we can conclude that $L_*^N \psi_{l, i}^- \leq 0$. Now we find that, $\phi_{l, 0}^- > 0, \; \phi_{l, N/2}^- > 0$ and $L^N_* \phi_{l, i}^- < 0$. Hence, by using Lemma 10 in ^[11], we can prove that $\phi_{l, i}^-\geq 0$. Therefore, $|W_{l, i}^-| \leq C \prod\limits_{j = 1}^{k_1}(1+\theta_2 h_j)^{-1}$.

Now consider the right layer barrier function $\psi_{r, i}^-.$ Applying the discrete operator defined in (3.11) over $\psi_{r, i}^-$, we obtain

Applying the central difference scheme in the place of $L^N_*$, we get

Using (4.1) for both the cases $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ and $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, we obtain

Similar arguments can be used to show that the upwind difference and mid-point difference satisfy

and

Hence, it is clear that, $\phi_{r, 0}^- > 0$, $\phi_{r, N/2}^- > 0$ and $L^N_* \phi_{r, i}^- < 0$. Therefore, applying Lemma 10 in ^[11], we obtain $\phi_{r, i}^- \geq 0$ and hence $|W_{r, i}^-|\leq C \prod\limits_{j = k_1+1}^{N/2}(1+\theta_1 h_j)^{-1}$. Similarly, we can prove the bounds of $W_{l, i}^+$ and $W_{r, i}^+$ as $i$ varies from $N/2+1$ to $N-1$. □

Now, we examine the truncation errors based on the three different operators. At the transition points, the mid-point scheme provides better accuracy in the convective terms compared to the central difference scheme on a non-uniform mesh. On a non-uniform mesh, we have

where $\overline{h}_i = (h_i + h_{i+1})/2$, and on a uniform mesh with step size h, we have

In the following lemma, we present the local error of the regular component.

Lemma 4.2. The regular component of the truncation error satisfies the following estimate:

Proof. For both the cases $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}\; \mbox{and}\; \sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d},$ when the mesh is uniform, i.e., $\tau_1 = \tau_2 = \tau_3 = \tau_4 = 1/8$, we have for $1 \leq i \leq N/2,$

Similarly, we obtain

When the mesh is non-uniform, we have

At the transition points, the upwind difference scheme is used if $2 \|b\|h_i \geq \varepsilon_c \alpha$ for $i = 2, 3, 5, 6$; otherwise, the mid-point difference scheme will be used. For both of these schemes, we obtain

Define the barrier function $\psi_i = CN^{-2} \zeta(x_i) \pm (V_i - v(x_i)),$ where

Then, $\psi_i$ satisfies

and $D^0 \psi_i \leq 0$, $D^+ \psi_i \leq 0$. Now, applying Lemma 3.1, we obtain

Therefore, we get the required result. □

The subsequent lemma provides the local error of the left and right singular components $W_{l, i}$ and $W_{r, i}$ associated with the boundary and interior layers.

Lemma 4.3. Let us assume (3.10). Then, the left and right singular components of the error satisfy the following estimates:

Proof. Firstly, consider the uniform mesh, i.e., $\tau_1 = \tau_2 = \tau_3 = \tau_4 = 1/8$, and for $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ for $i = 1, \cdots, N/2-1$, we have

Similarly, $|L^N_* (W_{l, i}^+ - w_l^+(x_i))| \leq C(N^{-1} \ln N)^2\; \forall \; N/2+1 \leq i \leq N-1.$

Again when $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, the above inequalities lead to

Similarly, $|L^N_* (W_{l, i}^+ - w_l^+)(x_i)| \leq C \varepsilon_c N^{-2} \ln^3 N.$

Now, we consider the error analysis on the non-uniform mesh. In the left boundary layer region $(0, \tau_1)$, the truncation error is

If $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, from (4.2) we obtain

Similarly, we can obtain the result $|L^N_* (W_{l, i}^+ - w_l^+(x_i))| \leq C (N^{-1} \ln N)^2$ on the left interior layer region $(d, d+\tau_3)$.

If $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$, we obtain from (4.2),

Now consider the barrier function on [$0, \tau_1$], see ^[11]

It is easy to check that the barrier function $\Psi_i$ is non-negative at the boundary points of [$0, \tau_1$] and $L_c^N \Psi_i < 0$. Hence, from Lemma 10 of ^[11], it follows $\Psi_i \geq 0$ for all [$0, \tau_1$]. Therefore

Similarly, we obtain $|W_{l, i}^+ - w_l^+(x_i)| \leq C N^{-2} \ln^3 N$ for $x_i \in (d, d+\tau_3)$.

When $x_i \in [\tau_1, d)$ for $\sqrt{\alpha} \varepsilon_d \leq \sqrt{\gamma \varepsilon_d}$ and $\sqrt{\alpha} \varepsilon_d \geq \sqrt{\gamma \varepsilon_d}$ we obtain the truncation error

Similarly, we have $|L^N_* (W_{l, i}^+ - w_l^+(x_i))| \leq C N^{-2}$ for $x_i \in [d+\tau_3, 1).$ Following the above procedures, we can prove the bound for the right singular components. □

The above estimates provide the truncation errors at the boundary and interior layer regions ($\Gamma^- \cup \Gamma^+$) except at the point of discontinuity. Estimating the error at the point of discontinuity is not straightforward. The following lemma gives an estimate for the local error at the point of discontinuity.

Lemma 4.4. At the point of discontinuity $x_{N/2} = d$, the error $e_d$ satisfies the following estimate:

Proof. Consider the case $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}.$ Then

For $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d},$ we have

where we have used Lemma 2.3 and a similar procedure adopted in ^[1]. □

Remark. When the sign of the discontinuous convection coefficient $a(x)$ is reversed, the above result at the point of discontinuity $x_i = d$ takes the form

5.   Error estimate

This section presents the main contribution of the paper. The following theorem provides the $\varepsilon_d$–$\varepsilon_c$ uniform higher-order error estimate of the computed solution.

Theorem 5.1. Let $u(x_i)$ be the solution of the continuous problem (1.1)-(1.2) and $U_i$ be the solution of the discrete problem (3.11) at $x = x_i$. Then, for sufficiently large $N$ satisfying the stability condition (3.10), we have

Proof. From the results of Lemmas 2.3, 4.2, 4.3 and using the procedure adopted in ^[11], it follows that

The presence of the discontinuity leads to the interior layers in a neighborhood of point $d$. We consider the following two cases to prove the error at the discontinuity point:

Case (i): $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$, define the discrete barrier function $\Phi_i$ to be the solution of the problem

where

A straight-forward application of the discrete minimum principle on the intervals $[0, d]$ and $[d, 1]$ leads to $0 \leq \Phi_i \leq 1.$ Also

Now, define the ancillary continuous functions $u_1(x)$ and $u_2(x)$ by

It is to be observed that the solutions of the above equations are

where

Define

Now

Hence, following the analysis given in ^[8] on the interval $[0, d]$ and $[d, 1]$, we obtain

and at the point of discontinuity, i.e., for $i = N/2$,

Now, define the mesh function

It can be clearly checked that $y_1(0), \; y_1(1) \geq 0, \; L^N_* y_1(x_i) \leq 0\; \forall \; x_i \in \overline{\Gamma}^N$ and $L^N_{t} y_1(x_{N/2}) \leq 0.$ Hence, Lemma 3.1 implies that $y_1(x_i) \geq 0\; \forall \; x_i \in \overline{\Gamma}^N$. Therefore

Case (ii): $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$. To show the error on the discontinuity point, we define the mesh function $y_2 (x_{i}) = \widetilde{W}_i \pm e_i$ on the mesh points in $(d-\tau_2, d+\tau_3)$ where,

Applying the central difference operator inside the domains $(d-\tau_2, d]$ and $[d, d+\tau_3)$, we obtain

Also,

Hence, we have

Now, $y_2 (d-\tau_2) \geq 0, \; y_2(d+\tau_3) \geq 0$, and $L^N_* y_2(x_i) \leq 0\; \mbox{for}\; x_i \in (d-\tau_2, d+\tau_3)$. Applying Lemma 3.1 to $y_2(x_i)$ in the above domain, we get

Hence, combining (5.1)–(5.3), we obtain the desired result. □

6.   Numerical examples

This section experimentally demonstrates the applicability of the hybrid scheme (3.11) and compares it with the existing upwind scheme (3.2) for two test problems. These problems have a jump discontinuity in the convective coefficient and source term. For these problems, the signs of the convection coefficients are different to show different boundary and interior layers. The numerical experiments are conducted on piecewise uniform Shishkin mesh, which changes with the various convection coefficients.

Example 6.1. Consider the two-parameter problem (1.1)-(1.2) with the following discontinuous convection coefficient and source term:

Example 6.2 Consider the singularly perturbed two-parameter BVP (1.1)-(1.2) with the following discontinuous convection coefficient and source term:

As the exact solutions of Examples 6.1 and 6.2 are unknown, we use the double mesh principle ^[5] to calculate the maximum pointwise error and the corresponding order of convergence of the numerical solution provided by the scheme (3.11). The error $E^N_{\varepsilon_d, \varepsilon_c}$ and corresponding order of convergence $\rho^N_{\varepsilon_d, \varepsilon_c}$ are computed as follows:

Here $U^N_i$ denotes the numerical solution obtained with $N$ number of mesh intervals, and $U^{2N}_i$ denotes the solution on $2N$ number of mesh intervals obtained by bisecting the previous original mesh. Similarly, we find the error $E^N$ and order of convergence $\rho^N$ of the existing upwind scheme (3.2) for a fixed value of $\varepsilon_c$ and various values of $\varepsilon_d,$ taken from the set $S = \{\varepsilon_d |\varepsilon_d = 10^{-2}, 10^{-4}, \cdots, 10^{-14}\}$:

Here, the numerical experiments are performed by choosing the constant values $\alpha$ = 1.25, $\beta$ = 1.0 and $\gamma$ = 0.8 for Example 6.1 and $\alpha$ = 2.0, $\beta$ = 1.0, and $\gamma$ = 0.5 for Example 6.2.

We present the maximum pointwise errors $E^N_{\varepsilon_d, \varepsilon_c}$ for Examples 6.1 and 6.2 in Tables 1 and 2, respectively. The corresponding orders of convergence $\rho^N_{\varepsilon_d, \varepsilon_c}$ for Examples 6.1 and 6.2 are shown in Tables 3 and 4, respectively. These tables show that the order of convergence is almost second-order $O(N^{-2} \ln^2 N)$ when $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ and $O(N^{-2} \ln^3 N)$ when $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$ for the above two examples. Table 5 presents the maximum error $E^N$ and order of convergence $\rho^N$ using the standard upwind scheme (3.2). This table shows that the existing scheme gives first-order parameter uniform convergence, while Tables 3 and 4 show almost second-order convergence using the hybrid difference scheme. Note that the errors are also lower in Tables 1 or 2 compared to Table 5.

Furthermore, we have plotted the numerical solution and the corresponding error for Example 6.1 in Figures 1 and 2. This shows that the interior layer (due to the discontinuity of the convection coefficient and source term) becomes sharper as $\varepsilon_c$ decreases from $10^{-2}$ to $10^{-4}$ for fixed $\varepsilon_d = 10^{-6}$ (i.e., the case $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ and $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$). Similar behavior can be observed for Example 6.2, whose solution and error graph are depicted in Figures 3 and 4. The Loglog plot of the maximum pointwise error for these problems in Figure 5 also demonstrates the uniform convergence of the numerical solution. From this figure we observe that the numerical (hybrid) method provides parameter-uniform convergence of $O(N^{-2} \ln^3 N)$, if $\sqrt{\alpha} \varepsilon_c \leq \sqrt{\gamma \varepsilon_d}$ and $O(N^{-2} \ln^2 N)$, if $\sqrt{\alpha} \varepsilon_c \geq \sqrt{\gamma \varepsilon_d}$ on Shishkin mesh.

7.   Conclusions

In this paper, an almost second-order uniformly convergent numerical solution is obtained for a two-parameter singularly perturbed problem where the convection coefficient and source term have a jump discontinuity at an interior point of the domain. The present hybrid difference scheme is a combination of upwind, mid-point, and central difference schemes at the interior points with a five-point scheme at the point of discontinuity so that it preserves the monotonicity property on the Shishkin mesh. Both theoretical and computational results based on this scheme produce almost second-order uniform convergence in the discrete maximum norm. Numerical experiments for the test problems validate the theoretical results. As an extension of this work, in the future, we aim to solve the proposed problem by considering a hybrid difference scheme using a Shishkin-Bakhvalov mesh.

Author contributions

All authors contributed equally to the writing of this article. All authors have accepted responsibility for the entire manuscript content and approved its submission.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in creating this article.

Acknowledgments

Vellore Institute of Technology, Vellore supported the work of MC under a SEED grant (Sanction Order No. SG20230081).

The authors would like to thank the editor and reviewers for their constructive comments.

Conflict of interest

Higinio Ramos is the(a) Guest Editor of special issue "Numerical Analysis of Differential Equations with Real-world Applications" for AIMS Mathematics. Higinio Ramos was not involved in the editorial review and the decision to publish this article.

The authors declare that there is no conflict of interest regarding the publication of this manuscript.