Cubic non-polynomial spline on piecewise mesh for singularly perturbed reaction differential equations with robin type boundary conditions

Abstract

Objective

The main purpose of this work is to present cubic non-polynomial spline approximation method for solving Robin-type singularly perturbed reaction–diffusion problems.

Results

The solution domain is first discretized using a piecewise mesh. The process begins by defining the cubic non-polynomial spline function and calculating its derivatives. These derivatives are then transformed into difference approximations, forming a linear system of algebraic equations in the form of a three-term recurrence relation, which is solved using an elimination algorithm. The stability and consistency of the method are analyzed, ensuring convergence. Numerical model examples are used to validate the proposed method, and the results are compared with those from other methods found in the literature. The maximum absolute error and the order of convergence for each example demonstrate the effectiveness and core contribution of the method.

Introduction

Mathematically modeled problems involve constructing a set of equations that consistently describe the characteristics or behavior of a physical system. Differential equations are often used to model various physical phenomena, where specific conditions are defined for the problem being studied [1]. A singularly perturbed differential equation arises when the highest order derivative is multiplied by a small positive parameter, called the perturbation parameter. These types of equations frequently appear in various fields of applied mathematics and engineering, including fluid mechanics, elasticity, quantum mechanics, chemical-reactor theory, aerodynamics, plasma dynamics, oceanography, meteorology, and the modeling of semiconductor devices [2,3,4,5].

The importance of these problems in real-world applications underscores the need for effective numerical methods to approximate their solutions. Numerical methods provide a practical approach when exact analytical solutions are either impossible or too complex to derive [2]. Among these, finite-difference methods are widely used for solving differential equations by approximating derivatives with finite differences. These methods transform differential equations into a system of algebraic equations that can be solved using iterative techniques.

Singularly perturbed problems are broadly categorized into reaction–diffusion and convection–diffusion types. As the perturbation parameter approaches zero, the order of the differential equation reduces by one or two [6,7,8,9]. These problems can also be classified based on boundary conditions, such as Dirichlet, Robin, or mixed conditions. In particular, second-order singularly perturbed reaction–diffusion problems with Robin boundary conditions exhibit boundary layers on both the left and right sides. This study focuses on the numerical solution of linear, second-order, singularly perturbed reaction–diffusion ordinary differential equations with Robin-type boundary conditions.

The challenge in solving singularly perturbed problems lies in the fact that the solution may vary rapidly in some parts of the domain and slowly in others, depending on the perturbation parameter [7, 8]. Standard numerical methods often fail to provide accurate results for these problems, especially as the perturbation parameter approaches zero. Therefore, it is crucial to develop numerical methods that are uniformly convergent and unaffected by the perturbation parameter. This study aims to develop such a method for Robin-type singularly perturbed reaction–diffusion problems.

Several numerical methods have been proposed in the literature for solving singularly perturbed reaction–diffusion problems [2, 5, 10,11,12,13,14], with many focusing on Dirichlet boundary conditions. However, more recent studies have addressed Robin-type boundary conditions [15,16,17,18], introducing higher-order finite difference methods for these problems. Also, the novel numerical methods, along with a detailed continuous solution analysis, are presented in ref. [15, 27, 28]. Despite these advances, there is still room for improvement in terms of the accuracy and convergence rates of the numerical solutions. Thus, the main objective of this work is to develop a higher-order, uniformly convergent numerical method for solving Robin-type singularly perturbed reaction–diffusion problems, along with a theoretical error analysis and numerical validation.

Description of the method

This paper deals with the Robin type singularly perturbed reaction diffusion problem of the form:

$$\left\{ \begin{aligned} Ly(x) \equiv & - \varepsilon y^{\prime\prime}(x) + b(x)y(x) = f(x),\,0 < x < 1, \\ B_{L} y(0) \equiv & a_{0} y(0) - a_{1} \sqrt \varepsilon y^{\prime}(0) = \eta_{1} , \\ B_{R} y(1) \equiv & d_{0} y(1) + d_{1} \sqrt \varepsilon y^{\prime}(1) = \eta_{2} , \\ \end{aligned} \right.$$

(1)

where $0 < {\varepsilon} < < 1$ is the singular perturbation parameter, and $a_{0} ,\,a_{1} ,\,d_{0} ,\,d_{1} ,\,\,\eta_{0} ,\,\eta_{1}$ with $a_{1} ,\,\,d_{1} \,\, \ne 0$ are given constants. The main purpose of this condition is to ensure the applicability of Robin-type boundary conditions. If this condition is not fulfilled, the problem shifts to a form with Dirichlet boundary conditions, which are more widely addressed by existing methods in the literature. The differential operators are denoted by $L,\,\,B_{L} ,\,\,B_{R}$.

The functions $b(x)$ and the source $f(x)$ function are assume to be sufficiently smooth functions such that $b(x) \ge b_{0} > 0,\,\,\,\forall x \in [0,\,\,1]$. Also, assume that the given constant values given in robin type boundary conditions satisfy $a_{j} ,\,\,d_{j} \,\, \ge 0$ and $a_{j} + \,\,d_{j} \,\, > 0,\,\,\,j = 0,1$.Then, the problem in Eq. (1) has a unique solution $y(x) \in C^{2} (0,\,\,1) \cap C^{1} [0,\,\,1]$ that typically exhibits boundary layers both at $x = 0$ and $x = 1$ when the perturbation parameter sufficiently small. The detailed proofs of the existence and uniqueness of the continuous solution for this defined problem is provided in ref. [15].

Piecewise mesh generation

We construct a piecewise mesh that contains more number of nodal points in the layer regions than non-layer region. The domain $[0,1],\,\,\,N \ge 4$ for $N$ is even multiple integer of 4, is divided into three subintervals, $[0,\tau ],\,\,[\tau ,1 - \tau ],\,\,[1 - \tau ,1]$ where the chosen transition parameter,

$$\tau = \min \left\{ {\frac{1}{4},\,\,\,2\sqrt \varepsilon \ln (N)} \right\},$$

(2)

denotes the width of the boundary layers. The domain $[0,1]^{N}$ is obtained by putting a non-uniform mesh with $\frac{N}{4}$ mesh elements in both the layer intervals and a uniform mesh with $\frac{N}{2}$ mesh elements in the outer layer region.

A general piecewise mesh $[0,1]^{N} = \left\{ {0 = x_{0} ,x_{1} ,x_{2} ,...,x_{N} = 1} \right\}$ with step size will be defined:

$$h_{i} = x_{i} - x_{i - 1} = \left\{ \begin{gathered} \frac{4\tau }{N},\,\,i = 1,2,...,\frac{N}{4} \hfill \\ \frac{2(1 - 2\tau )}{N},\,\,i = \frac{N}{4} + 1,...,\frac{3N}{4}, \hfill \\ \frac{4\tau }{N},\,\,i = \frac{3N}{4} + 1,...,N. \hfill \\ \end{gathered} \right.\,\,\,$$

(3)

Formulation of the method (cubic non-polynomial spline approximation)

Now, in order to develop the non-polynomial cubic spline approximation for the problems in Eq. (1), the interval $[0,\,\,1]$ is divided into N sub-intervals. Let $y(x_{i} )$ be the exact solution to the problem in Eq. (1), and $y_{i}$ be an approximation to $y(x_{i} )$, obtained by the segment $S(x)$ of the spline function passing through the points $(x_{i} ,\,\,y_{i} )$ and $(x_{i + 1} ,\,\,y_{i + 1} )$. Also, let us consider the non-polynomial cubic spline function $S(x)$ in subinterval $[x_{i} ,\,\,x_{i + 1} ],\,\,\,i = 1,\,\,...\,\,,\,N - 1$, of the form:

$$S(x) = c_{i} + d_{i} (x - x_{i} ) + u_{i} (e^{{\omega (x - x_{i} )}} - e^{{ - \omega (x - x_{i} )}} ) + v_{i} (e^{{\omega (x - x_{i} )}} + e^{{ - \omega (x - x_{i} )}} )$$

(4)

where, $c_{i} ,\,d_{i} ,\,u_{i} ,\,\,{\text{and}}\,\,v_{i}$ are unknown coefficients to be determined, and $\omega$ is a free parameter, which is used to raise the accuracy of the method.

To determine the unknown coefficient in Eq. (4), let us denote:

$$\begin{gathered} S(x_{i} ) = y_{i} ,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,S^{\prime}(x_{i} ) = m_{i} ,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,S^{\prime\prime}(x_{i} ) = M_{i} , \hfill \\ S(x_{i + 1} ) = y_{i + 1} ,\,\,\,\,\,\,\,\,\,S^{\prime}(x_{i + 1} ) = m_{i + 1} ,\,\,\,\,\,\,\,\,\,\,\,\,\,S^{\prime\prime}(x_{i + 1} ) = M_{i + 1} , \hfill \\ \end{gathered}$$

(5)

Differentiating Eq. (4) successively, we get:

$$S^{\prime}(x) = d_{i} + \omega u_{i} (e^{{\omega (x - x_{i} )}} + e^{{ - \omega (x - x_{i} )}} ) + \omega v_{i} (e^{{\omega (x - x_{i} )}} - e^{{ - \omega (x - x_{i} )}} )$$

(6)

$$S^{\prime\prime}(x) = \omega^{2} u_{i} (e^{{\omega (x - x_{i} )}} - e^{{ - \omega (x - x_{i} )}} ) + \omega^{2} v_{i} (e^{{\omega (x - x_{i} )}} + e^{{ - \omega (x - x_{i} )}} )$$

(7)

Using relations in Eqs. (5) and (7), we have:that results

$$\begin{gathered} S^{\prime\prime}(x_{i} ) = \omega^{2} u_{i} (e^{{\omega (x_{i} - x_{i} )}} - e^{{ - \omega (x_{i} - x_{i} )}} ) + \omega^{2} v_{i} (e^{{\omega (x_{i} - x_{i} )}} + e^{{ - \omega (x_{i} - x_{i} )}} ) = M_{i} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\, = \omega^{2} u_{i} (1 - 1) + \omega^{2} v_{i} (1 + 1) = M_{i} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\, = 2\omega^{2} v_{i} = M_{i} \hfill \\ \end{gathered}$$

$$v_{i} = \frac{{M_{i} }}{{2\omega^{2} }}.$$

(8)

Again, using the relation in Eqs. (5), (8) into Eq. (4) at the point $x_{i}$, yields:and it gives

$$\begin{gathered} S(x_{i} ) = c_{i} + d_{i} (x_{i} - x_{i} ) + u_{i} (e^{{\omega (x_{i} - x_{i} )}} - e^{{ - \omega (x_{i} - x_{i} )}} ) + v_{i} (e^{{\omega (x_{i} - x_{i} )}} + e^{{ - \omega (x_{i} - x_{i} )}} ) = y_{i} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\, = c_{i} + u_{i} (1 - 1) + v_{i} (1 + 1) = y_{i} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\, = \,c_{i} + 2v_{i} = y_{i} \hfill \\ \end{gathered}$$

$$c_{i} = y_{i} - \frac{{M_{i} }}{{\omega^{2} }}$$

(9)

Using the relation in Eqs. (5), (7) at the point $x_{i + 1}$ and letting $\theta_{2} = \omega h_{i + 1}$, we have:

$$\begin{gathered} S^{\prime\prime}(x_{i + 1} ) = \omega^{2} u_{i} (e^{{\omega (x_{i + 1} - x_{i} )}} - e^{{ - \omega (x_{i + 1} - x_{i} )}} ) + \omega^{2} v_{i} (e^{{\omega (x_{i + 1} - x_{i} )}} + e^{{ - \omega (x_{i + 1} - x_{i} )}} ) = M_{i + 1} \hfill \\ \hfill \\ \omega^{2} u_{i} (e^{{\omega h_{i + 1} }} - e^{{ - \omega h_{i + 1} }} ) + \omega^{2} v_{i} (e^{{\omega h_{i + 1} }} + e^{{ - \omega h_{i + 1} }} ) = M_{i + 1} \hfill \\ u_{i} (e^{{\theta_{2} }} - e^{{ - \theta_{2} }} ) + v_{i} (e^{{\theta_{2} }} + e^{{ - \theta_{2} }} ) = \frac{{M_{i + 1} }}{{\omega^{2} }} \hfill \\ u_{i} (e^{{\theta_{2} }} - e^{{ - \theta_{2} }} ) = \frac{{M_{i + 1} }}{{\omega^{2} }} - v_{i} (e^{{\theta_{2} }} + e^{{ - \theta_{2} }} ),\,\,\,\,\,\,\,\,\,\,\,\,\,But,\,\,v_{i} = \frac{{M_{i} }}{{2\omega^{2} }}\, \hfill \\ \end{gathered}$$

It gives

$$u_{i} = \frac{{M_{i + 1} }}{{\omega^{2} (e^{{\theta_{2} }} - e^{{ - \theta_{2} }} )}} - \frac{{M_{i} }}{{2\omega^{2} }}\frac{{(e^{{\theta_{2} }} + e^{{ - \theta_{2} }} )}}{{(e^{{\theta_{2} }} - e^{{ - \theta_{2} }} )}}$$

(10)

Also, using the relation in Eqs. (5), Eq. (4) at the point $x_{i + 1}$ and for $\theta_{2} = \omega h_{i + 1}$, we get:this results

$$\begin{gathered} S(x_{i + 1} ) = c_{i} + d_{i} (x_{i + 1} - x_{i} ) + u_{i} (e^{{\omega (x_{i + 1} - x_{i} )}} - e^{{ - \omega (x_{i + 1} - x_{i} )}} ) + v_{i} (e^{{\omega (x_{i + 1} - x_{i} )}} + e^{{ - \omega (x_{i + 1} - x_{i} )}} ) = y_{i + 1} \hfill \\ \hfill \\ y_{i + 1} = c_{i} + d_{i} h_{i + 1} + u_{i} (e^{{\omega h_{i + 1} }} - e^{{ - \omega h_{i + 1} }} ) + v_{i} (e^{{\omega h_{i + 1} }} + e^{{ - \omega h_{i + 1} }} ) \hfill \\ y_{i + 1} = c_{i} + d_{i} h_{i + 1} + u_{i} (e^{{\theta_{2} }} - e^{{ - \theta_{2} }} ) + v_{i} (e^{{\theta_{2} }} + e^{{ - \theta_{2} }} ) \hfill \\ \hfill \\ y_{i + 1} = y_{i} - \frac{{M_{i} }}{{\omega^{2} }} + d_{i} h_{i + 1} + \frac{{M_{i + 1} }}{{\omega^{2} }} - \frac{{M_{i} }}{{2\omega^{2} }}(e^{{\theta_{2} }} + e^{{ - \theta_{2} }} ) + \frac{{M_{i} }}{{2\omega^{2} }}(e^{{\theta_{2} }} + e^{{ - \theta_{2} }} ) \hfill \\ y_{i + 1} - y_{i} = d_{i} h_{i + 1} + \frac{{M_{i + 1} - M_{i} }}{{\omega^{2} }} \hfill \\ \hfill \\ \end{gathered}$$

$$d_{i} = \frac{{y_{i + 1} - y_{i} }}{{h_{i + 1} }} - \frac{{M_{i + 1} - M_{i} }}{{h_{i + 1} \omega^{2} }}$$

(11)

Using the continuity condition of the first derivative at $x_{i}$, $S^{\prime}_{\Delta - 1} (x_{i} ) = S^{\prime}_{\Delta } (x_{i} )$, we obtain:

$$d_{i - 1} + \omega u_{i - 1} (e^{{\theta_{1} }} + e^{{ - \theta_{1} }} ) - \omega v_{i - 1} (e^{{\theta_{1} }} - e^{{ - \theta_{1} }} ) = d_{i} + 2\omega u_{i}$$

(12)

Reducing indices of Eqs. (10) and (11) by one and substituting into Eq. (12), and after multiplying both sides by $\frac{2}{{h_{i} + h_{i + 1} }}$, and due to the parameters $\theta_{1} = \omega h_{i} ,\,\,\,{\text{and}}\,\,\,\theta_{2} = \omega h_{i + 1} ,$ we have the parameter $\omega = \frac{{\theta_{1} + \theta_{2} }}{{h_{i} + h_{i + 1} }}$, so that we get:

$$\begin{gathered} \frac{2}{{h_{i} + h_{i + 1} }}\left[ {\frac{{y_{i + 1} - y_{i} }}{{h_{i + 1} }} - \frac{{y_{i} - y_{i - 1} }}{{h_{i} }}} \right] = \alpha M_{i - 1} + 2\beta M_{i} \, + \alpha M_{i + 1} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered}$$

(13)

where the coefficients are denoted by:

$$\begin{gathered} \alpha = \frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{1}{{h_{i} \omega^{2} }} - \frac{{e^{{2\theta_{1} }} + e^{{ - 2\theta_{1} }} }}{{\omega (e^{{\theta_{1} }} - e^{{ - \theta_{1} }} )}}} \right) = \frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{1}{{h_{i + 1} \omega^{2} }} - \frac{2}{{\omega (e^{{\theta_{2} }} - e^{{ - \theta_{2} }} )}}} \right), \hfill \\ \beta = \frac{1}{{h_{i} + h_{i + 1} }}\left( {\frac{{(e^{{\theta_{1} }} + e^{{ - \theta_{1} }} )}}{{\omega (e^{{\theta_{1} }} - e^{{ - \theta_{1} }} )}} + \frac{{(e^{{\theta_{2} }} + e^{{ - \theta_{2} }} )}}{{\omega (e^{{\theta_{2} }} - e^{{ - \theta_{2} }} )}} - \frac{{h_{i} + h_{i + 1} }}{{h_{i} h_{i + 1} \omega^{2} }}} \right).\, \hfill \\ \end{gathered}$$

As $\theta_{1} \to 0,\,\,\,\theta_{2} \,\, \to 0$ in Eq. (13), we get the sum of the two constants, $\alpha + \beta = \frac{1}{2}$ and their values will be determined from local truncation error.

Considering the second order differential equation in Eq. (1) at the nodal point $x_{i}$ and using the relation in Eq. (5), corresponding to $M_{i} = y^{\prime\prime}_{i} = S^{\prime\prime}(x)$ which yields:

$$\begin{gathered} M_{i} = \frac{{b_{i} y_{i} - f_{i} }}{\varepsilon } \hfill \\ M_{i - 1} = \frac{{b_{i - 1} y_{i - 1} - f_{i - 1} }}{\varepsilon } \hfill \\ M_{i + 1} = \frac{{b_{i + 1} y_{i + 1} - f_{i + 1} }}{\varepsilon } \hfill \\ \end{gathered}$$

(14)

Thus, substituting Eq. (14) into Eq. (13) and also rearranging yields the finite difference scheme that re-written in three term recurrence relations:

$$E_{i} y_{i - 1} + F_{i} y_{i} + G_{i} y_{i + 1} = H_{i} ,\,\,\,\,\,\,\,\,i = 1,2,\,...\,,N - 1$$

(15)

where,

$$\begin{gathered} E_{i} = \frac{ - 2\varepsilon }{{h_{i} \left( {h_{i} + h_{i + 1} } \right)}} + \alpha b_{i - 1} \,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,F_{i} = \frac{2\varepsilon }{{h_{i} h_{i + 1} }} + 2\beta b_{i} , \hfill \\ G_{i} = \frac{ - 2\varepsilon }{{h_{i + 1} \left( {h_{i} + h_{i + 1} } \right)}} + \alpha b_{i + 1} ,\,\,\,\,\,\,\,\,\,\,\,\,\,H_{i} = \alpha f_{i - 1} + 2\beta f_{i} + \alpha f_{i + 1} . \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered}$$

Truncation error (To determine the values of $\alpha \,\,{\text{and}}\,\,\beta$)

For the formulated cubic spline method, the truncation error is obtained from Eq. (13) as:

$$\begin{gathered} T_{i} = \frac{2}{{h_{i} + h_{i + 1} }}\left[ {\frac{{y_{i + 1} - y_{i} }}{{h_{i + 1} }} - \frac{{y_{i} - y_{i - 1} }}{{h_{i} }}} \right] - (\alpha M_{i - 1} + 2\beta M_{i} \, + \alpha M_{i + 1} ) \hfill \\ \,\, \hfill \\ \end{gathered}$$

(16)

From the relation in Eq. (5) and by expanding Eq. (16) in Taylor’s series about $x_{i}$, we have the simplified truncation error:

$$\begin{gathered} T_{i} = \left( {1 - 2(\alpha + \beta )} \right)y^{\prime\prime}_{i} + \left( {h_{i + 1} - h_{i} } \right)\left\{ {\frac{1}{3} - \alpha } \right\}y^{\prime\prime\prime}_{i} + \left\{ {\frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{{h_{i + 1}^{3} + h_{i}^{3} }}{4!}} \right) - \frac{\alpha }{2}\left( {h_{i + 1}^{2} + h_{i}^{2} } \right)} \right\}y_{i}^{(4)} + \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\left\{ {\frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{{h_{i + 1}^{4} - h_{i}^{4} }}{5!}} \right) - \frac{\alpha }{3!}\left( {h_{i + 1}^{3} - h_{i}^{3} } \right)} \right\}y_{i}^{(5)} + \left\{ {\frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{{h_{i + 1}^{6} + h_{i}^{6} }}{6!}} \right) - \frac{\alpha }{4!}\left( {h_{i + 1}^{4} + h_{i}^{4} } \right)} \right\}y_{i}^{(6)} + ... \hfill \\ \end{gathered}$$

(17)

In order to obtain the higher order method, and for sufficiently cases of the mesh lengths, we choose the value of $\alpha$ and $\beta$ as:

$$\left\{ \begin{gathered} \alpha + \beta = \frac{1}{2},\, \hfill \\ \frac{1}{3} - \alpha = 0. \hfill \\ \end{gathered} \right.\,,\,\,\,\,\,\alpha = \frac{1}{3},\beta = \frac{1}{6},$$

$$\begin{aligned} \Rightarrow \frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{{h_{i + 1}^{3} + h_{i}^{3} }}{4!}} \right) - \frac{\alpha }{2}\left( {h_{i + 1}^{2} + h_{i}^{2} } \right) = & \frac{2}{{h_{i} + h_{i + 1} }}\left( {\frac{{h_{i + 1}^{3} + h_{i}^{3} }}{4!}} \right) - \frac{1}{6}\left( {h_{i + 1}^{2} + h_{i}^{2} } \right) \\ = & \frac{1}{6}\left\{ {\left( {\frac{{h_{i + 1}^{3} + h_{i}^{3} }}{{2(h_{i + 1} + h_{i} )}}} \right) - \left( {h_{i + 1}^{2} + h_{i}^{2} } \right)} \right\} \\ = & \frac{1}{6}\left\{ {\frac{{h_{i + 1}^{3} + h_{i}^{3} - 2(h_{i + 1}^{3} + h_{i} h_{i + 1}^{2} + h_{i + 1} h_{i}^{2} \, + h_{i}^{3} )}}{{2(h_{i + 1} + h_{i} )}}} \right\} \\ = & \frac{1}{6}\left\{ {\frac{{ - h_{i + 1}^{3} - 2h_{i} h_{i + 1}^{2} - 2h_{i + 1} h_{i}^{2} \, - h_{i}^{3} )}}{{2(h_{i + 1} + h_{i} )}}} \right\} \\ \le & \,(h_{i}^{2} + h_{i + 1}^{2} ) \\ \end{aligned}$$

and similarly, $\left\{ \begin{gathered} \alpha + \beta = \frac{1}{2},\, \hfill \\ \frac{1}{12} - \alpha = 0. \hfill \\ \end{gathered} \right.\,,\alpha = \frac{1}{12},\beta = \frac{5}{12},$

This yields the second order and fourth order convergent respectively.

Boundary conditions (end conditions)

For $i = 0,$ at $x_{0} = 0$, we have the discretized form of the first boundary condition, $B_{L} y(0)$ as

$$\,\,a_{0} y_{0} - a_{1} \sqrt \varepsilon y^{\prime}_{0} = \eta_{1}$$

(18)

In order to formulate the finite difference approximation to Eq. (18), adapting the Taylor series expansion,

$$y_{i + 1} = y_{i} + h_{i + 1} y^{\prime}_{i} + \frac{{h_{i + 1}^{2} }}{2}y^{\prime\prime}_{i} + \frac{{h_{i + 1}^{3} }}{3!}y^{\prime\prime\prime}_{i} \,\, + \frac{{h_{i + 1}^{4} }}{4!}y_{i}^{(4)} + \frac{{h_{i + 1}^{5} }}{5!}y_{i}^{(5)} + \frac{{h_{i + 1}^{6} }}{6!}y_{i}^{(6)} + ...$$

(19)

Considering Eq. (19) at the nodal point $x_{0}$, to approximation the first derivative in boundary condition, we obtain

$$y^{\prime}_{0} = \frac{{y_{1} - y_{0} }}{{h_{1} }} - \frac{{h_{1} }}{2}y^{\prime\prime}_{0} - \frac{{h_{1}^{2} }}{3!}y^{\prime\prime\prime}_{0} \,\, - \frac{{h_{1}^{3} }}{4!}y_{0}^{(4)} + TB_{L} ,$$

(20)

where the local truncation error for the left boundary condition is.

$TB_{L} = - \frac{{h_{1}^{4} }}{5!}y_{0}^{(5)} - \frac{{h_{1}^{5} }}{6!}y_{0}^{(6)} + ...$.

Differentiating the second order differential equation Eq. (1) twice with respect to the independent variable and considering it at the discretized points, we obtain

$$\begin{gathered} \,y^{\prime\prime}_{i} = \frac{{b_{i} y_{i} - f_{i} }}{\varepsilon } \hfill \\ y^{\prime\prime\prime}_{i} = \frac{{b_{i} y^{\prime}_{i} + b^{\prime}_{i} y_{i} - f^{\prime}_{i} }}{\varepsilon }, \hfill \\ \,y_{i}^{(4)} = \frac{{b_{i} y^{\prime\prime}_{i} + 2b^{\prime}_{i} y^{\prime}_{i} + b^{\prime\prime}_{i} y_{i} - f^{\prime\prime}_{i} }}{\varepsilon }. \hfill \\ \end{gathered}$$

(21)

Hence, from Eq. (21) the values of $y^{\prime\prime}_{0}$, $y^{\prime\prime\prime}_{0} \,\,,\,\,\,y_{0}^{(4)}$, obtained at the nodal point $x_{0}$ as:

$$\begin{gathered} y^{\prime\prime}_{0} = \frac{{b_{0} y_{0} - f_{0} }}{\varepsilon },\,\,\,\, \hfill \\ y^{\prime\prime\prime}_{0} = \frac{{b_{0} y^{\prime}_{0} + b^{\prime}_{0} y_{0} - f^{\prime}_{0} }}{\varepsilon }, \hfill \\ y_{0}^{(4)} = \frac{{b_{0} y^{\prime\prime}_{0} + 2b^{\prime}_{0} y^{\prime}_{0} + b^{\prime\prime}_{0} y_{0} - f^{\prime\prime}_{0} }}{\varepsilon }. \hfill \\ \end{gathered}$$

(22)

Substituting Eq. (22) into Eq. (20) and rearranging yields

$$y^{\prime}_{0} = \frac{{\frac{1}{{h_{1} }}y_{1} - \left\{ {\frac{{h_{1}^{2} b^{\prime}_{0} }}{3!\varepsilon } + \frac{{h_{1}^{3} b^{\prime\prime}_{0} }}{4!\varepsilon } + \frac{{b_{0} }}{\varepsilon }\left( {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right) + \frac{1}{{h_{1} }}} \right\}y_{0} + \frac{1}{\varepsilon }\left\{ {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right\}f_{0} + \frac{{h_{1}^{2} }}{3!\varepsilon }f^{\prime}_{0} + \frac{{h_{1}^{3} }}{4!\varepsilon }f^{\prime\prime}_{0} \,\,}}{{1 + \frac{{h_{1}^{2} b_{0} }}{3!\varepsilon } + \frac{{2h_{1}^{3} b^{\prime}_{0} }}{4!\varepsilon }}}.$$

(23)

Hence, for $i = 0,$ substituting Eq. (23) into Eq. (20), we get the scheme:

$$\begin{gathered} \left[ {a_{0} \left( {1 + \frac{{h_{1}^{2} b_{0} }}{3!\varepsilon } + \frac{{2h_{1}^{3} b^{\prime}_{0} }}{4!\varepsilon }} \right) + a_{1} \sqrt \varepsilon \left\{ {\frac{{h_{1}^{2} b^{\prime}_{0} }}{3!\varepsilon } + \frac{{h_{1}^{3} b^{\prime\prime}_{0} }}{4!\varepsilon } + \frac{{b_{0} }}{\varepsilon }\left( {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right) + \frac{1}{{h_{1} }}} \right\}} \right]y_{0} - \frac{{a_{1} \sqrt \varepsilon }}{{h_{1} }}y_{1} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = \eta_{1} \left( {1 + \frac{{h_{1}^{2} b_{0} }}{3!\varepsilon } + \frac{{2h_{1}^{3} b^{\prime}_{0} }}{4!\varepsilon }} \right) + a_{1} \sqrt \varepsilon \left\{ {\frac{1}{\varepsilon }\left\{ {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right\}f_{0} + \frac{{h_{1}^{2} }}{3!\varepsilon }f^{\prime}_{0} + \frac{{h_{1}^{3} }}{4!\varepsilon }f^{\prime\prime}_{0} \,} \right\}. \hfill \\ \end{gathered}$$

(24)

Similarly, for $i = N,$ at $x_{N}$, we have the discretized form of the second boundary condition, $B_{R} y(1)$ as

$$d_{0} y_{N} + d_{1} \sqrt \varepsilon y^{\prime}_{N} = \eta_{2} .$$

(25)

At any nodal point $x_{i}$, the Taylors’ series expansion written as

$$y_{i - 1} = y_{i} - h_{i} y^{\prime}_{i} + \frac{{h_{i}^{2} }}{2}y^{\prime\prime}_{i} - \frac{{h_{i}^{3} }}{3!}y^{\prime\prime\prime}_{i} \,\, + \frac{{h_{i}^{4} }}{4!}y_{i}^{(4)} - \frac{{h_{i}^{5} }}{5!}y_{i}^{(4)} + \frac{{h_{i}^{6} }}{6!}y_{i}^{(6)} + ...$$

(26)

Then, adapting this expansion for the nodal point $x_{N}$, to approximation the first derivative in the second boundary condition, we get:

$$y^{\prime}_{N} = \frac{{y_{N} - y_{N - 1} }}{{h_{N} }} + \frac{{h_{N} }}{2}y^{\prime\prime}_{N} - \frac{{h_{N}^{2} }}{3!}y^{\prime\prime\prime}_{N} \,\, + \frac{{h_{N}^{3} }}{4!}y_{N}^{(4)} + TB_{R} ,$$

(27)

where the local truncation error for the second boundary condition is

$TB_{R} = - \frac{{h_{N}^{4} }}{5!}y_{N}^{(5)} + \frac{{h_{N}^{5} }}{6!}y_{N}^{(6)} + ...$.

Similar to Eq. (22), from Eq. (21) at the nodal point $x_{N}$, considering the values of $y^{\prime\prime}_{N}$, $y^{\prime\prime\prime}_{N} \,\,,\,\,\,y_{N}^{(4)}$, and also substituting it into Eq. (27), we get:

$$y^{\prime}_{N} = \frac{{\frac{{ - y_{N - 1} }}{{h_{N} }} + \left\{ {\frac{1}{{h_{N} }} - \frac{{h_{N}^{2} b^{\prime}_{N} }}{3!\varepsilon } + \frac{{h_{N}^{3} b^{\prime\prime}_{N} }}{4!\varepsilon } + \frac{{h_{N} b_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N}^{2} }}{{4!\varepsilon^{2} }}} \right\}y_{N} - \left( {\frac{{h_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N} }}{{4!\varepsilon^{2} }}} \right)f_{N} + \frac{{h_{N}^{2} f^{\prime}_{N} }}{3!\varepsilon }\,\, - \frac{{h_{N}^{3} f^{\prime\prime}_{N} }}{4!\varepsilon }}}{{1 + \frac{{h_{N}^{2} b_{N} }}{3!\varepsilon } - \frac{{2h_{N}^{3} b^{\prime}_{N} }}{4!\varepsilon }}}.$$

(28)

Now, substituting this value of $y^{\prime}_{N}$ in Eq. (28) into Eq. (25) gives the scheme

$$\begin{gathered} \frac{{ - d_{1} \sqrt \varepsilon }}{{h_{N} }}y_{N - 1} + \left[ {d_{0} \left( {1 + \frac{{h_{N}^{2} b_{N} }}{3!\varepsilon } - \frac{{2h_{N}^{3} b^{\prime}_{N} }}{4!\varepsilon }} \right) + d_{1} \sqrt \varepsilon \left( {\frac{1}{{h_{N} }} - \frac{{h_{N}^{2} b^{\prime}_{N} }}{3!\varepsilon } + \frac{{h_{N}^{3} b^{\prime\prime}_{N} }}{4!\varepsilon } + \frac{{h_{N} b_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N}^{2} }}{{4!\varepsilon^{2} }}} \right)} \right]y_{N} \hfill \\ \,\,\,\,\,\, \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = \eta_{2} \left( {1 + \frac{{h_{N}^{2} b_{N} }}{3!\varepsilon } - \frac{{2h_{N}^{3} b^{\prime}_{N} }}{4!\varepsilon }} \right) + d_{1} \sqrt \varepsilon \left( {\frac{{h_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N} }}{{4!\varepsilon^{2} }}} \right)f_{N} - d_{1} \sqrt \varepsilon \frac{{h_{N}^{2} f^{\prime}_{N} }}{3!\varepsilon }\,\, + d_{1} \sqrt \varepsilon \frac{{h_{N}^{3} f^{\prime\prime}_{N} }}{4!\varepsilon }. \hfill \\ \end{gathered}$$

(29)

To conclude formulation of the method, the three equations in Eqs. (15), (24) and (29) constitute the system of $(N + 1) \times (N + 1)$ linear algebraic equations, that give the approximations $y_{0} ,\,\,y_{1} ,\,\,...\,\,,\,y_{N}$ of the solution $y(x)$ at the nodal points,$x_{0} ,\,\,x_{1} ,\,\,...\,\,,\,x_{N}$, with piecewise mesh length $h_{i}$. This can be re-written in matrix form

$$MY = R,$$

(30)

where the three matrices are defined

$M = (m_{ij} )_{(N + 1) \times (N + 1)} = \left[ {\begin{array}{*{20}c} {B_{0} } & {C_{0} } & 0 & \cdots & 0 \\ {A_{1} } & {B_{1} } & {C_{1} } & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & {A_{N - 1} } & {B_{N - 1} } & {C_{N - 1} } \\ 0 & \cdots & 0 & {A_{N} } & {B_{N} } \\ \end{array} } \right],\,\,\,\,\,\,\,\,\,\,Y = \left[ {\begin{array}{*{20}c} {y_{0} } \\ {y_{1} } \\ \vdots \\ {y_{N - 1} } \\ {y_{N} } \\ \end{array} } \right],\,\,\,\,\,\,\,\,\,\,R = \left[ {\begin{array}{*{20}c} {H_{0} } \\ {H_{1} } \\ \vdots \\ {H_{N - 1} } \\ {H_{N} } \\ \end{array} } \right]$,

with the specific entries are denoting:

$$\begin{gathered} {\text{For}}\,\,\,i = 0,\,\,\,\,B_{0} = a_{0} \left( {1 + \frac{{h_{1}^{2} b_{0} }}{3!\varepsilon } + \frac{{2h_{1}^{3} b^{\prime}_{0} }}{4!\varepsilon }} \right) + a_{1} \sqrt \varepsilon \left( {\frac{{h_{1}^{2} b^{\prime}_{0} }}{3!\varepsilon } + \frac{{h_{1}^{3} b^{\prime\prime}_{0} }}{4!\varepsilon } + \frac{{b_{0} }}{\varepsilon }\left( {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right) + \frac{1}{{h_{1} }}} \right), \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,C_{0} = - \frac{{a_{1} \sqrt \varepsilon }}{{h_{1} }}, \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,H_{0} \, = \eta_{1} \left( {1 + \frac{{h_{1}^{2} b_{0} }}{3!\varepsilon } + \frac{{2h_{1}^{3} b^{\prime}_{0} }}{4!\varepsilon }} \right) + a_{1} \sqrt \varepsilon \left\{ {\frac{1}{\varepsilon }\left( {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right)f_{0} + \frac{{h_{1}^{2} }}{3!\varepsilon }f^{\prime}_{0} + \frac{{h_{1}^{3} }}{4!\varepsilon }f^{\prime\prime}_{0} \,} \right\}, \hfill \\ \end{gathered}$$

$${\text{For}}\,\,i = 1,\,\,...\,\,,\,\,N - 1:\,\,\,A_{i} \, = \,E_{i} ,\,\,\,\,B_{i} \, = \,F_{i} ,\,\,\,\,C_{i} \, = \,G_{i} ,\,\,\,\,\,\,{\text{and}}\,\,\,\,\,H_{i} \, = \,\alpha f_{i - 1} + 2\beta f_{i} + \alpha f_{i + 1} .$$

$$\begin{gathered} {\text{For}}\,\,\,i = N:\,\,\,\,A_{N} = \frac{{ - d_{1} \sqrt \varepsilon }}{{h_{N} }}, \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,B_{N} = d_{0} \left( {1 + \frac{{h_{N}^{2} b_{N} }}{3!\varepsilon } - \frac{{2h_{N}^{3} b^{\prime}_{N} }}{4!\varepsilon }} \right) + d_{1} \sqrt \varepsilon \left( {\frac{1}{{h_{N} }} - \frac{{h_{N}^{2} b^{\prime}_{N} }}{3!\varepsilon } + \frac{{h_{N}^{3} b^{\prime\prime}_{N} }}{4!\varepsilon } + \frac{{h_{N} b_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N}^{2} }}{{4!\varepsilon^{2} }}} \right), \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,H_{N} \, = \eta_{2} \left( {1 + \frac{{h_{N}^{2} b_{N} }}{3!\varepsilon } - \frac{{2h_{N}^{3} b^{\prime}_{N} }}{4!\varepsilon }} \right) + d_{1} \sqrt \varepsilon \left( {\frac{{h_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N} }}{{4!\varepsilon^{2} }}} \right)f_{N} - d_{1} \sqrt \varepsilon \frac{{h_{N}^{2} f^{\prime}_{N} }}{3!\varepsilon }\,\, + d_{1} \sqrt \varepsilon \frac{{h_{N}^{3} f^{\prime\prime}_{N} }}{4!\varepsilon }. \hfill \\ \end{gathered}$$

Algorithm and its stability

The obtained finite difference approximation yields an $(N + 1) \times (N + 1)$ algebraic tri-diagonal system that can be solved by Gaussian elimination method. Hence, in using this method, we have the algorithm:

$${\text{Forward}}\,\,{\text{Elimination}}\,\,\,\, \to \,\,\left\{ \begin{gathered} W_{0} = \frac{{C_{0} }}{{B_{0} }}, \hfill \\ T_{0} = \frac{{R_{0} }}{{B_{0} }}, \hfill \\ W_{i} = \frac{{C_{i} }}{{B_{i} - A_{i} W_{i - 1} }},\,\,\,\,\,i = 1,\,2,\,\,\,...\,\,,\,N + 1, \hfill \\ T_{i} = \frac{{R_{i} - A_{i} T_{i - 1} }}{{B_{i} - A_{i} W_{i - 1} }},\,\,\,\,\,i = 1,\,2,\,\,\,...\,\,,\,N + 1. \hfill \\ \end{gathered} \right.$$

$${\text{Backward}}\,\,{\text{Elimination}}\,\,\,\, \to \,\,\left\{ \begin{gathered} {\text{Setting}}\,\,\,C_{N} = 0, \hfill \\ y_{N + 1} = T_{N + 1} , \hfill \\ y_{i} = T_{i} - W_{i} y_{i + 1} ,\,\,\,\,\,i = \,N,\,\,N - 1,\,\,...\,\,\,2,\,\,1. \hfill \\ \end{gathered} \right.$$

Theorem 1: (Stability)

Let $M$ be a coefficient matrix of the tri-diagonal system, in Eq. (30). Then, the matrix $M$ is an irreducible and diagonally dominant matrix and hence the scheme is stable [26].

Proof

To estimate the stability of the formulated method, we show that the coefficient matrix of the method satisfy the diagonal dominance. First, to check this, consider that $b(x) \ge b_{0} > 0,\,\,\,\forall x \in [0,\,\,1]$ and for $i = 1,\,\,...\,\,,\,\,N - 1$, from Eq. (16), we have:

$$\begin{gathered} \,B_{i} \,\, \ge \,\,A_{i} + C_{i} ,\,\,\,\,(F_{i} \,\, \ge \,\,E_{i} + G_{i} ),\,\,\,\,\,\,\,i = 1,2,\,\,\,...\,\,,\,N - 1, \hfill \\ \frac{2\varepsilon }{{h_{i} h_{i + 1} }} + 2\beta b_{i} \,\,\, \ge \frac{ - 2\varepsilon }{{h_{i} \left( {h_{i} + h_{i + 1} } \right)}} + \alpha b_{i - 1} \, + \frac{ - 2\varepsilon }{{h_{i + 1} \left( {h_{i} + h_{i + 1} } \right)}} + \alpha b_{i + 1} , \hfill \\ \frac{2\varepsilon }{{h_{i} h_{i + 1} }} + 2\beta b_{0} \,\,\, \ge \frac{ - 2\varepsilon }{{h_{i} \left( {h_{i} + h_{i + 1} } \right)}} + \alpha b_{0} \, + \frac{ - 2\varepsilon }{{h_{i + 1} \left( {h_{i} + h_{i + 1} } \right)}} + \alpha b_{0} , \hfill \\ \frac{2\varepsilon }{{h_{i} h_{i + 1} }} + \frac{{5b_{0} }}{6}\,\,\, > \frac{ - 2\varepsilon }{{h_{i} h_{i + 1} }} + \frac{{b_{0} }}{12}. \hfill \\ \end{gathered}$$

Second, we consider the Robin boundary conditions at the two extremities $i = 0$ and $i = N.$

$$\begin{gathered} {\text{For}}\,\,\,i = 0,\,\,\,\,B_{0} \, > \,\,\,C_{0} , \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,B_{0} = a_{0} \left( {1 + \frac{{h_{1}^{2} b_{0} }}{3!\varepsilon } + \frac{{2h_{1}^{3} b^{\prime}_{0} }}{4!\varepsilon }} \right) + a_{1} \sqrt \varepsilon \left( {\frac{{h_{1}^{2} b^{\prime}_{0} }}{3!\varepsilon } + \frac{{h_{1}^{3} b^{\prime\prime}_{0} }}{4!\varepsilon } + \frac{{b_{0} }}{\varepsilon }\left( {\frac{{h_{1}^{3} b_{0} }}{4!\varepsilon } + \frac{{h_{1} }}{2}} \right) + \frac{1}{{h_{1} }}} \right),\,\,\,\,\,\,C_{0} = - \frac{{a_{1} \sqrt \varepsilon }}{{h_{1} }}. \hfill \\ \end{gathered}$$

$$\begin{gathered} {\text{For}}\,\,\,i = N:\,\,\,\,\,\,\,B_{N} \,\, > \,\,\,\,A_{N} , \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,A_{N} = \frac{{ - d_{1} \sqrt \varepsilon }}{{h_{N} }},\,\,\,B_{N} = d_{0} \left( {1 + \frac{{h_{N}^{2} b_{N} }}{3!\varepsilon } - \frac{{2h_{N}^{3} b^{\prime}_{N} }}{4!\varepsilon }} \right) + d_{1} \sqrt \varepsilon \left( {\frac{1}{{h_{N} }} - \frac{{h_{N}^{2} b^{\prime}_{N} }}{3!\varepsilon } + \frac{{h_{N}^{3} b^{\prime\prime}_{N} }}{4!\varepsilon } + \frac{{h_{N} b_{N} }}{2\varepsilon } + \frac{{h_{N}^{3} b_{N}^{2} }}{{4!\varepsilon^{2} }}} \right). \hfill \\ \end{gathered}$$

Thus, this implies that for each row of $M$, the sum of the two off-diagonal elements is less than or equal to the modulus of the diagonal element. Therefore, $M$ is diagonally dominant. Hence, $M$ is irreducible matrix. Henceforth, the Thomas Algorithm (Solution of general Tri-diagonal system) or Gaussian elimination method for solving the formulated scheme in Eq. (30) is stable for the described method.

Numerical examples and results

In order to test the validity of the proposed method and to demonstrate their convergence computationally, we have taken different model examples. The maximum absolute error is evaluated by:

$$E_{\varepsilon }^{N} = \mathop {\max }\limits_{{x_{i} \in [0,1]^{N} }} |y\left( {x_{i} } \right) - y_{i} |,\,\,$$

where $y\left( {x_{i} } \right)$ and $y_{i}$ respectively, denotes the exact and approximate solutions. Further, the order of convergence calculated by the formula:

$$R = \frac{{\log (E_{\varepsilon }^{N} ) - \log (E_{\varepsilon }^{2N} )}}{\log (2)}.\,\,\,\,$$

Example 1

Consider the singularly perturbed problem [15]

$$\left\{ \begin{gathered} - \varepsilon y^{\prime\prime}(x) + (1 + x(1 - x))y(x) = f(x),\,\,\,\,\,\,0 < x < 1, \hfill \\ \,\,\,\,\,y(0) - \sqrt \varepsilon y^{\prime}(0) = - 1 - \sqrt \varepsilon + e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} , \hfill \\ \,\,\,\,\,y(1) + \sqrt \varepsilon y^{\prime}(1) = - 1 - \sqrt \varepsilon + e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} . \hfill \\ \end{gathered} \right.$$

The source function is given by

$$f(x) = 1 + x(1 - x) + [2\sqrt \varepsilon - x(1 - x)^{2} ]e^{{\frac{ - x}{{\sqrt \varepsilon }}}} + [2\sqrt \varepsilon - x^{2} (1 - x)]e^{{\frac{ - (1 - x)}{{\sqrt \varepsilon }}}} .$$

The exact solution to the problem that satisfies the above boundary conditions is given by.

$y(x) = 1 + (x - 1)e^{{\frac{ - x}{{\sqrt \varepsilon }}}} - xe^{{\frac{ - (1 - x)}{{\sqrt \varepsilon }}}}$.

Example 2

Consider the singularly perturbed problem [16]

$$\left\{ \begin{gathered} - \varepsilon y^{\prime\prime}(x) + \frac{4(1 + \sqrt \varepsilon (1 + x))}{{(1 + x)^{4} }}y(x) = f(x),\,\,\,\,\,\,0 < x < 1, \hfill \\ \,\,\,\,\,y(0) - \sqrt \varepsilon y^{\prime}(0) = 2 + \frac{6}{{1 - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}, \hfill \\ \,\,\,\,\,y(1) + \sqrt \varepsilon y^{\prime}(1) = \frac{{ - 2 - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}{{2(1 - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} )}}, \hfill \\ \end{gathered} \right.$$

With the source function is $f(x) = \frac{ - \,4}{{(1 + x)^{4} }}$$\left[ {(1 + } \right.\sqrt \varepsilon (1 + x) + 4\pi^{2} \varepsilon )\cos (\frac{4\pi x}{{1 + x}})$$- 2\pi \varepsilon (1 + x)\sin (\frac{4\pi x}{{1 + x}}) + 3(1 + \sqrt \varepsilon (1 + x))\,\frac{{e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}{{1 - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}].$

The mixed boundary conditions are chosen in a way, such that the exact solution for this problem becomes

$y(x) = - \cos (\frac{4\pi x}{{1 + x}}) + 3\frac{{e^{{\frac{ - 2x}{{(1 + x)\sqrt \varepsilon }}}} - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}{{1 - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}$.

Example 3

Consider the singularly perturbed problem [15]

$$\left\{ \begin{gathered} - \varepsilon y^{\prime\prime}(x) + y(x) = - \cos^{2} (\pi x) = 2\varepsilon \pi^{2} \cos (2\pi x),\,\,\,\,\,\,0 < x < 1, \hfill \\ \,\,\,\,\,y(0) - 2\sqrt \varepsilon y^{\prime}(0) = 2 - \frac{{4e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}{{1 + e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}, \hfill \\ \,\,\,\,\,y(1) + 3\sqrt \varepsilon y^{\prime}(1) = \frac{{3(1 - e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} )}}{{1 + e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }}. \hfill \\ \end{gathered} \right.$$

The exact solution that satisfies the given boundary conditions is given by

$$y(x) = \frac{{e^{{\frac{ - (1 - x)}{{\sqrt \varepsilon }}}} + e^{{\frac{ - x}{{\sqrt \varepsilon }}}} }}{{1 + e^{{\frac{ - 1}{{\sqrt \varepsilon }}}} }} - \cos^{2} (\pi x)$$

Example 4

Consider singularly perturbed problem [16]

$$\left\{ \begin{gathered} - \varepsilon y^{\prime\prime}(x) + y(x) = - \cos^{2} (\pi x) - 2\varepsilon \pi^{2} \cos (2\pi x),\,\,\,\,\,\,x \in (0,1) \hfill \\ \,\,\,\,\,y(0) - 2\sqrt \varepsilon y^{\prime}(0) = 1, \hfill \\ \,\,\,\,\,y(1) + 3\sqrt \varepsilon y^{\prime}(1) = 0. \hfill \\ \end{gathered} \right.$$

As the exact solution for the Example 4 is not available, so the accuracy of its numerical solution will be computed using double mesh principle (Figs. 1, 2 and Tables 1, 2, 3, 4, 5, 6, 7).

Fig. 1

Effects of perturbation parameter on the solution for Example 1 when, $N = 64$

Fig. 2

Numerical solution versus exact solution for Examples 2 and 3 respectively, when $N = 64,\,\,\varepsilon = 10^{ - \,4}$

Table 1 Computed maximum absolute errors for Example 1 by the proposed method on uniform mesh when $\alpha = \frac{1}{12},\,\,\beta = \frac{5}{12}$

Table 2 Computed maximum absolute errors for Example 1 by the proposed method on piecewise mesh

Table 3 Computed maximum absolute errors for Example 2 by the proposed method

Table 4 Computed maximum order of convergence for Example 2 by the proposed method

Table 5 Comparison of maximum absolute errors and orders of convergence for Example 2, on the set for $\varepsilon \in S = \left\{ {2^{ - 2} ,\,\,...\,\,,\,\,2^{ - \,40} } \right\}$

Table 6 Comparison of maximum absolute errors and orders of convergence for Example 3

Table 7 Comparison of maximum absolute errors and orders of convergence for Example 4, on the set for $\varepsilon \in S = \left\{ {2^{ - 2} ,\,\,...\,\,,\,\,2^{ - \,40} } \right\}$

Conclusion

In conclusion, this paper presents a novel family of non-polynomial cubic spline approximation methods for solving singularly perturbed reaction–diffusion problems with Robin-type boundary conditions. By utilizing piecewise mesh lengths, this approach effectively addresses the challenge posed by boundary layers, resulting in a numerical scheme of either second- or fourth-order convergence. Through the transformation of boundary conditions into a fourth-order finite difference scheme and the construction of a tri-diagonal system, the method achieves stability and consistency, ensuring convergence. The diagonally dominant nature of the coefficient matrix allows for efficient solution using Gaussian elimination, highlighting the method’s computational effectiveness.

The numerical results further emphasize the advantages of this method, particularly in handling the small positive perturbation parameter that compromises accuracy on a uniform mesh. When applied to a piecewise mesh, finer discretization’s in the boundary layer and coarser sizes in the outer region significantly improve accuracy, as demonstrated by numerical illustrations. The decrease in maximum absolute error with increasing mesh points substantiates the method's convergence. The proposed cubic spline approach outperforms existing methods in terms of accuracy and convergence, closely approximating the exact solution with fewer errors. Overall, this method offers a highly accurate and efficient solution for Robin-type singularly perturbed problems, aligning well with theoretical predictions.

Limitations

Various approaches can be applied for the finite difference approximation of differential equations and boundary conditions. However, these differing methods may introduce inaccuracies when applied to the specific problem at hand.