Fitted mesh methods based on non-polynomial splines for singularly perturbed boundary value problems with mixed shifts

1.   Introduction

Delay differential equations model problems in several domains, including biosciences, material science, and medicine^[1,2]. Differential-difference equations are differential equations in which the system's evolution depends on both the system's historical context and its present state. A differential equation that consists of at least one shift term and whose highest-order derivative is multiplied by a small perturbation parameter is known as a singularly perturbed differential-difference equation. Singularly perturbed differential difference equations (SPDDEs) generally lead to solutions exhibiting boundary layers, and as the perturbation parameter goes to zero, the smoothness of the solution deteriorates.

The initial developments on the asymptotic analysis of singularly perturbed differential-difference equations have emerged from the articles by Lange and Miura ^[3,4]. The numerical explorations in this field can be found in ^[5], where the authors Kadalbajoo and Sharma have presented a $\epsilon$-uniform numerical scheme comprising of a standard upwind finite difference operator on a fitted piecewise uniform mesh for a class of boundary value problems of singularly perturbed differential-difference equations with small shifts. In ^[6], the authors Kadalbajoo et al. developed a fitted operator and a fitted mesh finite difference method for a class of singularly perturbed difference-difference equations, the fitted mesh method being $\epsilon$-uniform convergent of second order. Kadalbajoo and Ramesh ^[7] proposed a hybrid numerical method on Shishkin mesh for solving a singularly perturbed delay differential equation, wherein the solutions were compared with those obtained by using a simple upwind scheme and a midpoint upwind scheme. Sirisha et al. ^[8] presented a mixed finite difference method for singularly perturbed differential-difference equations with mixed shifts via domain decomposition on a constant mesh. Woldaregay and Duressa ^[9] presented a numerical scheme for singularly perturbed differential-difference equations with mixed small shifts by using the exponentially fitted operator finite difference method, and in ^[10] they applied an exponentially finite difference method to a singularly perturbed boundary value problem. Ranjan and Prasad ^[11] used an exponentially fitted three-term finite difference technique to approximate the solution of a singularly perturbed differential equation with small shifts. Kumar and Kadalbajoo ^[12] constructed a piecewise uniform mesh for solving singularly perturbed differential-difference equations with small shifts.

Jain^[13] introduced the spline function approximation and shown that they are the consistency relations for the fundamental equations in discrete mechanics. Kadalbajoo and Bawa ^[14] proposed a variable mesh difference scheme for singularly perturbed boundary value problems using splines, and the method is shown to be quadratically convergent. Kadalbajoo and Patidar ^[15,16] applied spline techniques such as spline in compression and spline in tension to singularly perturbed two-point boundary value problems. Aziz and Khan ^[17] studied a spline method for second-order singularly perturbed boundary value problems, and the convergence of the method is shown to be dependent on the choice of the parameters. Mohanty and Jha ^[18] applied the variable mesh method using spline in compression for singularly perturbed two-point singular boundary value problems. Mohanty and Arora ^[19] applied tension spline on a non-uniform mesh for singularly perturbed two-point singular boundary value problems with significant first derivatives. Chakravarthy et al.^[20,21] presented the numerical solution using spline in compression and spline in tension on a uniform mesh for second-order singularly perturbed delay differential equations. Ravi Kanth and Murali ^[22] presented a numerical technique for solving singularly perturbed nonlinear delay differential equations by the method of Spline in compression. The quasilinearization technique is applied in converting the nonlinear equation into a sequence of linear equations.

The aforementioned publications illustrate the implementation of spline methods for various types of singularly perturbed boundary value problems and motivate us in exploring the feasibility of constructing non-polynomial spline methods for solving singularly perturbed boundary value problems with mixed shifts. In the present paper, we constructed numerical methods using spline in compression, spline in tension and adaptive spline on a fitted mesh, and a comparative study is performed on the results. This method ensures a consistent level of accuracy regardless of the perturbation parameter and achieves reliable convergence.

The content of this paper is organized as follows: In Section 2, we introduce the problem under consideration and some properties of the solution. In Section 3, we discuss the mesh construction strategy for the problem. Section 4 is devoted to the proposed methods for the problem. In Section 5, we discuss the convergence of the proposed methods. In Section 6, we present the numerical results for test problems, and finally the conclusions follow in Section 7.

2.   Statement of the problem

Consider the general boundary value problem (BVP) for SPDDE containing both types of shifts:

for $x\in (0, 1), 0 < \epsilon\ll1$, subject to the interval and boundary conditions

where $0 < \epsilon\ll1$ is the perturbation parameter, $a(x), \alpha(x), \omega(x), \beta(x), f(x)$, $\phi(x)$, and $\chi(x)$ are smooth functions, and $\delta$ is the delay (negative shift) parameter and $\eta$ is the advance (positive shift) parameter. As $\delta, \eta < \epsilon$, for $a(x)-\delta \alpha(x)+\eta \beta(x) > 0$ and $\alpha(x)+\omega(x)+\beta(x) < 0$ $\forall$ $x\in [0, 1]$, the solution exhibits a boundary layer near $x = 0$, while for $a(x)-\delta \alpha(x)+\eta \beta(x) < 0$ and $\alpha(x)+\omega(x)+\beta(x) < 0$, the solution exhibits a boundary layer near $x = 1$. Here we assume that $\alpha(x)\leq M_1$, $\beta(x) \leq M_2$, and $\alpha(x)+\omega(x)+\beta(x) \le - \mathscr{M} < 0$. The function $y(x)$ being continuous in $[0, 1]$ and differentiable in $(0, 1)$, satsifying (2.1) and (2.2), provides a smooth solution for (2.1) and (2.2).

Since the solution of $y(x)$ of (2.1) and (2.2) is sufficiently differentiable, we expand the terms $y(x-\delta)$ and $y(x-\eta)$ using Taylor series to obtain:

On substituting Eq (2.3) in Eq (2.1), the modified form of the Eq (2.1) is

subject to the conditions

where $\Upsilon(x) \approx y(x)$, $\mu(x) = \epsilon+\alpha(x)\frac{\delta^2}{2}+\beta(x)\frac{\eta^2}{2}$, $p(x) = a(x)-\delta \alpha(x)+\eta \beta(x)$, $q(x) = \alpha(x)+\omega(x)+\beta(x)$ and $r(x) = f(x)$.

2.1.   Some properties of the solution

We show that the operator $\mathscr{L}$ follows the minimum principle for the continuous problem Eq (2.4):

Lemma 1. Let $\Upsilon(x)$ be a smooth function, with $\Upsilon(0)\geq 0, \Upsilon(1)\geq 0$, then for $x\in [0, 1]$, $\Upsilon(x) \geq 0$, whenever $\mathscr{L} (\Upsilon(x))\leq 0$ for $x\in (0, 1)$.

Proof. For proof of this lemma, the reader can refer to ^[5]. □

The bound for the solution of the continuous problem (2.4) is given in the following lemma:

Lemma 2. Let $\Upsilon(x)$ be the solution of (2.4) and (2.5), then, $\Vert \Upsilon \Vert \leq \mathscr{M}^{-1} \Vert r \Vert + \max (\vert \phi_0 \vert, \vert \chi_1 \vert)$, $\Vert. \Vert$ being the $l_\infty$ norm $\Vert \Upsilon \Vert = \max_{s\in[0, 1]} \vert \Upsilon(s) \vert$.

Proof. For proof of this lemma, the reader can refer to ^[5]. □

Lemma 1 guarantees the uniqueness of the solutions of (2.4) and (2.5), and the existence of the solution is guaranteed as the given problem is linear. Also, the boundedness to the solution of the problem is implied by Lemma 2. Also, the bounds for the solutions of (2.4) and (2.5) and their derivatives are given in the following lemma.

Lemma 3. Considering $\Upsilon(x)$ to be the solution of (2.4) and (2.5), we have $\vert \vert \Upsilon^{(k)} \vert \vert \leq 2^kC\left(2\epsilon+ \delta^2 M_1+ \eta^2 M_2\right)^{-k}$ for $k = 1, 2, 3$.

Proof. For proof of this lemma, the reader can refer to ^[5]. □

Theorem 1. Let $\Upsilon(x)$ be the solution of (2.4) and (2.5), and let $\Upsilon(x) = \Upsilon_r(x)+\Upsilon_s(x)$, where the regular component $\Upsilon_r(x)$ satisfies

and the singular component $\Upsilon_s(x)$satisfies

where $0\leq k\leq3$.

Proof. For proof of this theorem, the reader can refer to ^[5]. □

3.   Mesh construction

In this section, we discuss the mesh generation for the numerical solution of the singularly perturbed BVP (2.4) and (2.5).

The case when the boundary layer occurs at the left end of the domain $D = [0, 1]$, it is divided into two subdomains $D_1$ and $D_2$ such that $D = D_1\cup D_2 = [0, \tau] \cup [\tau, 1]$, where $\tau$ is the transition parameter that is closer to $x = 0$, and is defined by

where $N$ is the number of mesh points in the domain $D = [0, 1]$ and $\tau_0\geq \frac{1}{\vert \mathscr{M} \vert}$. It is clear that $\tau = \frac{1}{2}$, the mesh is uniform; otherwise, the mesh condenses near the left boundary. It is assumed that $N = 2^m$, where $m\geq2$ is an integer, which guarantees that there is at least one point in the boundary layer region. So, we consider equal number of mesh points in each subdomain and uniform partition over each subdomain with mesh points $x_i$, as defined by

where $h_1 = \frac{2\tau}{N}$ and $h_2 = \frac{2(1-\tau)}{N}$ on the domains $D_1$ and $D_2$ respectively.

Similarly, in the case when the boundary layer occurs at the right end of the solution domain $D$, we divide into subdomains $D^*_1$ and $D^*_2$ such that $D = D^*_1\cup D^*_2 = [0, 1-\tau] \cup [1-\tau, 1]$, where $\tau$ is so-called the transition parameter and is located near the point $x = 1$. We consider equal number of grid points in each subdomain and uniform partition over each subdomain with grid points $x_i$, as defined by

Now we show that $\mathscr{L}$ satisfies the discrete minimum principle:

Lemma 4. If the mesh function $\Upsilon(x_i)$ satisfying $\Upsilon(x_0)\geq 0, \Upsilon(x_N)\geq 0$, then $\Upsilon(x_i)\geq 0$, $0\leq x_i \leq 1$, for $\mathscr{L} (\Upsilon(x_i))\leq 0$, $0 < x_i < 1$.

Proof. Let $0\leq \bar z_k \leq 1$ be such that $\Upsilon(\bar z_k) = \min_{x\in[0, 1]} \Upsilon(x_i)$, and assuming that $\Upsilon(\bar z_k) < 0$, clearly $\bar z_k\notin \{0, 1 \}$. Hence $\Upsilon^\prime(\bar z_k) = 0$ and $\Upsilon^{\prime \prime}(\bar z_k) \geq 0.$

Now we have $\mathscr{L}(\Upsilon(\bar z_k)) = \mu(\bar z_k)\Upsilon^{\prime \prime}(\bar z_k)+p(x_i) \Upsilon^\prime(\bar z_k)) + q(x_i)\Upsilon(\bar z_k) > 0,$ which is a contradiction to our assumption that $\Upsilon(\bar z_k) < 0$. Therefore, $\Upsilon(\bar z_k) \geq 0$ and hence $\Upsilon(x_k) \geq 0 \; \forall x_i \in[0, 1]$. □

Lemma 5. Let $\Upsilon_i$ be any mesh function such that $\Upsilon_0 = \Upsilon_N = 0$. Then, for all $0\leq i\leq N$, $\Vert \Upsilon_j \Vert \leq \mathscr{M}^{-1} \max \limits_{1\leq j \leq N-1} \vert \mathscr{L}(\Upsilon_j) \vert$.

Proof. Let us introduce two mesh functions $\widehat{\upsilon_i}^\pm$ defined by

and for $1\leq i \leq N-1$

Therefore, by the discrete minimum principle, we have $\widehat{\upsilon}_i^\pm \geq 0 \ \forall i, 0\leq i\leq N$, which gives the required estimate. □

Lemma 6.

for each $i$.

Proof.

The above inequality is true for each $j$. Now we multiply these inequalities for $j = i+1, ..., N$, and we obtain

Hence the result. □

Lemma 7. For $i = 0, 1, ..., N$, we set

then for $i = 0, 1, ..., N-1$, we have

Proof. It is easy to verify that

Now

from which the result follows. □

Lemma 8. There exists a constant $C$ such that

for $N/2\leq i \leq N$.

Proof. suppose $N/2\leq i \leq N$. By ^[23]

It is easy to verify that $N^{16(i-1/N)N^{-1}\ln N /(1+4N^{-1}\ln N)}$ is bounded for any $N\geq2$ from which the result follows. □

4.   Numerical methods

In this section, we present non-polynomial spline methods for solving the boundary value problems (2.4) and (2.5).

4.1.   Spline in compression

The spline in compression $S_{\Delta}(x)$ satisfies in $[x_{i-1}, x_{i}]$ the differential equation

where $S_{\Delta}(x_{i}) = \Upsilon_{i}, \psi > 0, h_i = x_i-x_{i-1}$.

Solving (4.1) as a second-order differential equation, we obtain

Applying the interpolating conditions at $x_{i-1}$ and $x_{i}$; $S_{\Delta}(x_{i-1}) = \Upsilon_{i-1}$, $S_{\Delta}(x_{i}) = \Upsilon_{i}, S^{\prime\prime}_{\Delta}(x_{i}) = M_{i}$, and setting $\lambda_i = \sqrt \psi h_i$, we obtain the interpolating constants $A$ and $B$ and hence

Differentiating (4.3) and taking $x\to x_i$, we obtain

Considering the interval $(x_i, x_{i+1})$ and similarly we obtain

Equating the left and right hand dertivatives at $x_{i}$, we obtain

This leads to a tridiagonal system

where $\overline\lambda = \frac{1}{\lambda_{i}^2}\left[\frac{\lambda_{i}}{\sin \lambda_{i}}-1\right]$ and ${\overline{\overline\lambda}} = \frac{1}{\lambda_{i}^2}[1-\lambda_{i} \cot \lambda_{i}]$.

The consistency relation for (4.5) may be expressed as $\frac{\lambda_i}{2} = \tan(\frac{ \lambda_i}{2})$, whose smallest positive root is $\lambda_i \approx 8.986818916$, which leads to the equation ${\overline\lambda}+{\overline{\overline\lambda}} = \frac{1}{2}$.

To obtain an approximation for $\Upsilon_i^\prime$ and $\Upsilon_i^{\prime\prime}$, we use the Taylor series approximation for $\Upsilon$ about $x_i$ as:

Solving (4.6) and (4.7) for $\Upsilon_i^\prime$ and $\Upsilon_i^{\prime\prime}$, we will obtain

Using the above approximations (4.8) and (4.9) in $\Upsilon_{i+1}^\prime = \Upsilon_i^\prime+h_{i+1}\Upsilon_i^{\prime\prime}$ and $\Upsilon_{i-1}^\prime = \Upsilon_i^\prime-h_i \Upsilon_i^{\prime\prime}$, we obtain

We write Eq (2.4) as

Now we rewrite Eq (4.5) as

Substituting (4.12) in (4.13) and using the approximations (4.8), (4.10), and (4.11), we obtain the following tridiagonal scheme:

where

Using the Thomas algorithm, we can solve the above tri-diagonal scheme (4.14) subject to the boundary conditions (2.5).

4.2.   Spline in tension

The spline in tension $S_{\Delta}(x)$ in $[x_{i-1}, x_{i}]$ satisfies the differential equation

where $\psi > 0$ is a tension factor, $S_{\Delta}(x_{i}) = \Upsilon_{i}$, $S^\prime_{\Delta}(x_{i}) = m_{i}$, $S^{\prime\prime}_{\Delta}(x_{i}) = M_{i}$, $h_i = x_i-x_{i-1}$.

Solving (4.15) as a second-order differential equation, we obtain

Applying the interpolating conditions at $x_{i-1}$ and $x_{i}$ and setting $\Lambda_i = \sqrt \psi h_i$, we obtain the interpolating constants $A$ and $B$, and hence

Differentiating (4.17) and taking $x\to x_i$, we obtain

Considering the interval $(x_i, x_{i+1})$ and similarly, we obtain

Equating the left and right hand dertivatives at $x_{i}$, we have

This leads to a tridiagonal system

where $\Lambda_{1} = \frac{1}{\Lambda_i^2}\left[1-\frac{\Lambda_i}{\sinh \Lambda_i}\right]$ and $\Lambda_{2} = \frac{1}{\Lambda_i^2}[\Lambda_i \coth \Lambda_i-1]$.

We rewrite the Eq (4.19) as

Substituting (4.12) in (4.20) and using the approximations (4.8), (4.10), and (4.11), we obtain the following tridiagonal linear system:

where

Using the Thomas algorithm, we can solve the above tri-diagonal scheme (4.21) subject to the boundary conditions (2.5).

4.3.   Adaptive spline

The function $S_{\Delta}(x)$, which we call adaptive spline, satisfies the following differential equation:

Solving (4.22) and using interpolatory constraints $S_{\Delta}(x_{i-1}) = \Upsilon_{i-1}$, $S_{\Delta}(x_{i}) = \Upsilon_{i}$, we obtain

where $\nu_i = \frac{\psi h_i}{2\Theta}$, $z_i = \frac{x-x_{i-1}}{h_i}$ and $\Theta, \psi$ are constants.

The function $S_{\Delta}(x)$ on the interval $[x_i, x_{i+1}]$ is obtained by replacing $i$ with $i+1$ (4.23).

Applying the conditions of continuity to the first or second derivative of $S_{\Delta}(x)$ at $x_i$, we obtain the following relationship:

which simplifies to the following form of tridiagonal system:

Some additional relations for the adaptive spline are listed as follows:

(i) $m_{i-1} = -h_i(A_1M_{i-1}+A_2M_i)+\frac{\Upsilon_{i}-\Upsilon_{i-1}}{h_i}$

(ⅱ) $m_i = h_i(A_3M_{i-1}+A_4M_i)+\frac{\Upsilon_{i}-\Upsilon_{i-1}}{h_i}$

(ⅱ) $M_{i-1} = \frac{2\nu_i}{\varsigma_i h_i}\left[-(A_4m_{i-1}+A_2m_i)+B_1(\frac{\Upsilon_{i}-\Upsilon_{i-1}} {h_i})\right]$

(ⅳ) $M_i = \frac{2\nu_i}{\varsigma_i h_i}\left[(A_3m_{i-1}+A_4m_i)+B_2(\frac{\Upsilon_{i}-\Upsilon_{i-1}}{h_i})\right]$

where $\varsigma_i = \frac{\nu_i coth \nu_i-1}{2\nu_i}$,

$A_1 = \frac{1}{4}(1+2\varsigma_i)+\frac{\varsigma_i}{2\nu_i}$, $A_2 = \frac{1}{4}(1-2\varsigma_i)-\frac{\varsigma_i}{2\nu_i}$, $A_3 = \frac{1}{4}(1+2\varsigma_i)-\frac{\varsigma_i}{2\nu_i}$ and

$A_4 = \frac{1}{4}(1-2\varsigma_i)+\frac{\varsigma_i}{2\nu_i}$, $B_1 = \frac{1}{2}(1-2\varsigma_i)$, $B_2 = \frac{-1}{2}(1+2\varsigma_i)$.

In the limiting case, when $\nu_i \to 0$, we have

$\varsigma_i = 0$, $\frac{\varsigma_i}{\nu_i} = \frac{1}{6}$, $A_1 = \frac{1}{3}$, $A_2 = \frac{1}{6}$, $A_3 = \frac{1}{6}$, $A_4 = \frac{1}{3}$, $B_1 = \frac{1}{2}$, $B_2 = \frac{-1}{2}$

and the spline function (4.23) reduces to the cubic spline.

By introducing the parameter $\mu_i$, we rewrite Eq (4.27) as

Substituting (4.12) in (4.28) and using the approximations (4.8), (4.10), and (4.11), we obtain the tridiagonal linear system of the form:

where

Using the Thomas algorithm, we can solve the above tri-diagonal scheme (4.29) subject to the boundary conditions (2.5).

5.   Convergence analysis

Here we perform the convergence analysis for the scheme described in Section 4.1.

Writing the tri-diagonal system Eq (4.14) in matrix-vector form, we obtain

in which $A = [m_{i, j}], 1\leq i, j \leq N-1$, is a tri-diagonal matrix of order $N-1$, with

and $C = (d_i)$ is a column vector with $d_i = h_i^2(\overline\lambda r_{i-1}+2\overline{\overline\lambda}r_{i}+\overline\lambda r_{i+1})$ with $i = 1, 2, \cdots, N-1$ with $T(h_i) = O(h_i^3)$.

We also have

where $(\bar \Upsilon) = (\bar{\Upsilon_0}, \bar{\Upsilon_1}, \cdots, \bar{\Upsilon_N})^T$ and $T(h_i) = (T_0(h_i), T_1(h_i), \cdots, T_N(h_i))^T$ denote the actual solution and the local truncation error, respectively.

From Eqs (5.1) and (5.2), we obtain

Thus the error equation is

where $E = \bar{\Upsilon}-\Upsilon = (e_o, e_1, e_2, \cdots, e_N)^T$. Let $\vert p(x)\vert \leq c_1$, $\vert q(x)\vert \leq c_2$ and $[m_{i, j}]$ is the $(i, j)^{th}$ element of the matrix $A$. Then we have

For sufficiently small $h_i$, we have

Hence the matrix is irreducible ^[24].

Let the $i^{th}$ row elements' sum of matrix $A$ be $S_i$, then we have

Let $c_{1^*} = \min\vert p(x)\vert$, $c^*_1 = \max\vert p(x)\vert$, $c_{2^*} = \min\vert q(x)\vert$, $c^*_2 = \max\vert q(x)\vert$ and $h = \max\limits_{i = 1}^{N-1} \{h_i, h_{i+1}\}$ so that $0 < c_{1^*}\leq c_1 \leq c^*_1, 0 < c_{2^*}\leq c_2 \leq c^*_2$.

Then for a given $h$, the matrix $A$ is irreducible and monotone (^[24,25]).

From (5.3), we have

At the end points $i = 0$ and $N$, the above inequality holds, and for $1\leq i \leq N-1$, we have

where $\hat\zeta\in(x_{i-1}, x_i)$.

Since the mesh is piecewise uniform with step difference $h$, from (5.5) and (5.6), we obtain

Also, we have from ^[26]

as $0 < h < 1$.

Using (5.7) and (5.8) in (5.5), we obtain

For a left layer problem, let the fine mesh points for the inside layer region be $x_1, \cdots, x_{N/2}$, and the coarse mesh points in the outer region be $x_{N/2+1}, \cdots, x_{N-1}$. Further $\Upsilon_r(x)$ and $\Upsilon_s(x)$ are the regular and the singular components of the numerical solution. Then, using (5.2) along with (3.1) and $h_1$, $h_2$ into (5.9) and Theorem 1, we obtain

Similarily,

The above results (5.10) and (5.11) can be concluded as

Theorem 2. $a(x), \alpha(x), \omega(x), \beta(x), f(x)$, $\phi(x)$, and $\chi(x)$ be sufficiently smooth functions so that $\Upsilon(x)\in C^3 [0, 1]$. Let $\Upsilon_i, \ i = 0(1)N$ be the approximate solution of (2.4), obtained using the fitted mesh finite difference method (4.14) with the conditions (2.5). Then, there is a constant $\mathcal{M}$ independent of $\epsilon$ and the mesh size such that

6.   Numerical results

To check the efficiency of the methods described in Sections 4.1–4.3, four test problems of SPDDEs are solved, of which two problems are of left end boundary layer type and the other two are right layer problems.

The double mesh principle is used for finding the maximum absolute errors, which is given by the formula:

and the numerical rate of convergence for the considered problems is calculated by the following formula:

The numerical techniques outlined in Sections 4.1–4.3 are applied to the test problems, and the maximum absolute errors and the numerical rate of convergence are evaluated. The numerical results are tabulated for a spectrum of values of $\delta$ and $\eta$, smaller than $\epsilon$ for all the test problems. The findings are displayed in Tables 1–8. Also, the $\epsilon$-uniform maximum absolute errors $E_N$ for various values of the mesh parameter $N$ and for $\epsilon \in \{2^0, 2^{-1}, \cdots 2^{-20}\}$ is compared for each method described in Section 4 in Table 9. The numerical work illustrates the efficiency of the methods and is also consistent with those in literature. Graphs illustrating the influence of the shift parameters on the solution of the problem are depicted in Figures 1–8. The relationship between the error $E_N$ and the number of mesh points $N$ for the considered examples is plotted in Figures 9–12. These plots illustrate the efficiency of the methods presented in Sections 4.1–4.3.

Example 1. $\epsilon y^{\prime\prime}(x)+ y^{\prime}(x)-2y(x-\delta)+y(x)-y(x+\eta) = -1$, $y(x) = 1;-\delta \leq x\leq 0, y(x) = 1, 1\leq x \leq 1+\eta$.

Example 2. $\epsilon y^{\prime\prime}(x)+2.5y^{\prime}(x)-2e^{x}y(x-\delta)-y(x)-xy(x+\eta) = 1$, $y(x) = 1;-\delta \leq x\leq 0, y(x) = 1, 1\leq x \leq 1+\eta$.

Example 3. $\epsilon y^{\prime\prime}(x)- y^{\prime}(x)-2y(x-\delta)+y(x)-2y(x+\eta) = 0$, $y(x) = 1;-\delta \leq x\leq 0, y(x) = -1, 1\leq x \leq 1+\eta$.

Example 4. $\epsilon y^{\prime\prime}(x)-(1+e^{x^2})y^{\prime}(x)-xy(x-\delta)+x^2y(x)-(1-e^{-x})y(x+\eta) = 1$, $y(x) = 1;-\delta \leq x\leq 0, y(x) = -1, 1\leq x \leq 1+\eta$.

7.   Discussion and conclusions

In this paper, we proposed fitted mesh numerical methods for solving singularly perturbed boundary value problems of second-order ordinary differential equations with mixed shifts. The shifts that are smaller than the perturbation parameter are approximated using Taylor series and non-polynomial splines, namely, the spline in compression, the spline in tension, and the adaptive spline are applied to the Shishkin mesh. The methods presented are analyzed for convergence and shown to be first-order convergent. Numerical computations are carried out on two test problems that exhibit layer behavior on the left of the underlying interval and two right-layer problems. The maximum absolute errors and rates of convergence are tabulated, which show the first-order uniform rate of convergence. Graphs are plotted for the test problems for different values of the perturbation and shift parameters. From the figures, the effect of the shifts on the boundary layer behavior of the solution of the problems can be observed. As the shifts increase in magnitude, the thickness of the layer decreases for the left-layer problems, while for the right-layer problems, it increases. The methods presented in this paper have been found to be almost equally efficient in achieving $\epsilon$-uniform convergence and the numerical rate of convergence. Hence, it can be concluded that the presented methods provide considerable advantage for solving singularly perturbed linear second-order boundary value problems with mixed shifts.

Author contributions

Both the authors contributed equally to this work.

Acknowledgments

The authors wish to thank the National Board for Higher Mathematics, Department of Atomic Energy, Government of India, for their financial support under the project No. 02011/8/2021 NBHM(R.P)/R&D Ⅱ/7224, dated 24.06.2021.

Conflict of interest

The authors declare that they have no conflict of interest, relevant to the content of this article.