Parameter-Robust Numerical Scheme for Computing Singularly Perturbed Differential Equations with Mixed Large Shifts

1. Introduction

Singularly perturbed functional differential equations (SP-FDEs) are a type of differential equation that combines the features of functional differential equations (FDEs) and singularly perturbed differential equations. They involve both mixed-shift arguments and small perturbation parameters. Functional differential equations are a class of differential equations where the unknown function depends not only on the current values of its variables but also on past values. FDEs include time delays or memory effects in the equation. The solutions of FDEs often require the history of the function over a certain time interval [1]. Singularly perturbed differential equations are also a class of differential equations where a small parameter $\epsilon$ appears as a coefficient of the term with largest order of derivative. This small parameter characterizes the relative importance of the different terms in the equation [2]. The presence of the mixed-shift arguments introduces memory effects, where the solution at a given time depends not only on the current state but also on its history. The small parameter $\epsilon$ determines the relative importance of the different terms in the equation. As $\epsilon$ approaches zero, the differential equation becomes singularly perturbed (i.e. it becomes increasingly stiff) with rapid variations in some regions of the solution.

SP-FDEs are arise in many applications across different scientific disciplines. They offer a powerful framework for modeling systems with memory effects and perturbed dynamics. For instance, in control theory, the observed in modeling systems with delays, memory effects, and perturbed dynamics [3]. In population dynamics, it arise in modeling the interactions, delays, and effects of environmental changes [4]. In epidemiology, studying the spread of infectious diseases, incorporating time delays in transmission rates are modeled by SP-FDEs [5]. In biology and bio-engineering, it is frequently observed in modeling biological systems with memory effects, such as gene regulation, enzyme kinetics, or drug delivery [6].

Solving singularly perturbed functional differential equations poses several challenges due to the combined complexities of functional dependencies, delay effects, and small perturbation parameters. These problems often exhibit multi-scale behavior, where the solution varies at different rates in different regions of the domain. This can lead to sharp gradients, boundary or interior layers, and slow-fast dynamics. Capturing and resolving these variations accurately is a challenge. Obtaining analytical solutions for SP-FDEs is exceedingly difficult due to the combination of functional dependencies, delays, and perturbation parameters. Applying analytical techniques to general SP-FDEs are limited, and solving equations with specific forms or simplifying assumptions may be required. On the other hand due to the small parameter $\epsilon$, standard numerical methods can struggle with stability and convergence issues. The presence of delays and singular perturbations requires carefully designed numerical schemes that can handle the specific characteristics of SP-FDEs. Addressing these difficulties often requires specialized knowledge of numerical analysis, stability theory, delay differential equations, and the characteristics of functional dependencies.

To overcome these challenges, various researchers have developed different numerical schemes using a combination of numerical methods, adaptive algorithms, and careful analysis of the problem’s behavior and solution properties. The choice of the most appropriate approach depends on the specific SP-FDE being solved and its unique characteristics. For instance, in [7], a singularly perturbed problems with shift is solved by formulating a fitted exponential spline method. The authors accomplished derivation of the method by including a fitted factor in the constructed numerical scheme which is based on a special type of grids, and the fitted factor is evaluated applying the theory of singular perturbation. Woldaregay [8] solved singularly perturbed differential equations with delays by developing three different schemes. The author used nonstandard mid-point upwind finite difference method on uniform grids, a standard mid-point upwind finite difference method on Shiskin’s meshes and a nonstandard mid-point upwind finite difference method on Shiskin’s meshes. Numerical analysis and results of the methods show that the third technique is more accurately convergent method.

In the work of [9] a non-polynomial cubic spline technique is devised to solve singularly perturbed delay problems having two layer regions and oscillatory characteristics in the solution. Lalu and Phaneendra [10] developed a difference scheme for nonlinear singularly perturbed differential equation using the trigonometric spline technique and obtained a uniformly convergent method. In [11], singularly perturbed differential equations involving mixed large delays of convection-diffusion type are solved by developing a fitted mesh numerical method. The authors used finite and hybrid difference schemes to treat the problem on piece-wise uniform Shishkin mesh. Elango and Unyong [12] solved singularly perturbed SP-FDEs involving mixed delays. They treated the problem by constructing fitted mesh scheme on Shishkin-type meshes and Bakhvalov-Shishkin type meshes. In the work of [13], a fitted operator scheme is constructed by decomposing the domain into inner and outer parts for singularly perturbed ordinary differential equations and obtained a second order convergent method. Wondimu et al. [14] treated singularly perturbed problems with delay and non-local boundary condition by developing an exponentially fitted finite difference scheme and obtained that the method is uniformly convergent.

Though various forms of perturbation problems have been treated in literature, very few attentions are given for singularly perturbed differential equations with mixed large shifts. In this paper, motivated by the fitted mesh methods developed in [11,12], we inspired to develop a fitted operator difference scheme for the reaction-diffusion type singularly perturbed ordinary differential equations with mixed large shifts. To tackle the influence of the shift arguments, we have chosen a special uniform mesh in which the shift arguments lie on the nodal points. On such meshes, a numerical scheme is proposed using the Numerov’s finite difference approach.The convergence analysis of the scheme shows that the method is parameter-robust numerical scheme.

The remaining sections of this work are arranged as follows: In Section 2, the continuous problem and some preliminary properties of the continuous solution are presented. The numerical formulation and analysis are presented in Section 3. Numerical examples are provided in Section 4 to confirm the theoretical findings in practice. In Section 5 some concluding remarks and future plans end the study in this paper.

Notation: In this work, $\mathtt{C}$ is taken as a generic number not depending on the mesh number and the perturbation parameter. The norm $\parallel .\parallel$ is defined as $\parallel .\parallel ={max}_{\overline{\Phi }}|.|$ for any function in $\overline{\Phi }$.

2. Continuous Problem

Consider the singularly perturbed functional differential equation given by

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }\omega \left(x\right)=& -\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}+\mu \left(x\right)\omega \left(x\right)+\nu \left(x\right)\omega (x-1)+\gamma \left(x\right)\omega (x+1)=\psi \left(x\right),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}\omega \left(x\right)=p\left(x\right),x\in [-1,0],\hfill \\ \omega \left(x\right)=q\left(x\right),x\in [3,4],\hfill \end{array}\right.\hfill \end{array}$$

(2)

where $0<\epsilon \ll 1$ and ${\mathcal{L}}_{\epsilon }$ is a differential operator. The coefficient functions $\mu \left(x\right)$, $\nu \left(x\right)$, $\gamma \left(x\right)$ and the source function $\psi \left(x\right)$ are assumed to be continuous on $[0,3]$, and the history functions $p\left(x\right)$ and $q\left(x\right)$ are continuous on $[-1,0]$ and $[3,4]$, respectively. To avoid oscillations in the solution, we assumed that the coefficient functions satisfy

$$\mu \left(x\right)+\nu \left(x\right)+\gamma \left(x\right)\ge 2\alpha >0,\mu \left(x\right)>0,\nu \left(x\right)<0and\gamma \left(x\right)<0,\forall x\in [0,3].$$

(3)

For simplicity, we used the notations ${\Phi }^{-}=(0,1)$, ${\Phi }^0=(1,2)$, ${\Phi }^{+}=(2,3)$, $\Phi =(0,3)$, $\overline{\Phi }=[0,3]$ and ${\Phi }^{*}={\Phi }^{-}\cup {\Phi }^0\cup {\Phi }^{+}$. Thus, the above assumptions hold that for $x\in \overline{\Phi }$, $\omega \in \mathtt{M}={\mathtt{C}}^0\left(\overline{\Phi }\right)\cap {\mathtt{C}}^1(\Phi )\cap {\mathtt{C}}^2({\Phi }^{-}\cup {\Phi }^0\cup {\Phi }^{+})$.

Considering the interval conditions in (2), we rewrite (1) as

$${\mathcal{L}}_{\epsilon }\omega \left(x\right)=\Re \left(x\right),$$

(4)

where

$$\begin{array}{cc}\hfill & {\mathcal{L}}_{\epsilon }\omega \left(x\right)=\left\{\begin{array}{c}-\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}+\mu \left(x\right)\omega \left(x\right)+\gamma \left(x\right)\omega (x+1),x\in {\Phi }^{-},\hfill \\ -\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}+\mu \left(x\right)\omega \left(x\right)+\nu \left(x\right)\omega (x-1)+\gamma \left(x\right)\omega (x+1),x\in {\Phi }^0,\hfill \\ -\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}+\mu \left(x\right)\omega \left(x\right)+\nu \left(x\right)\omega (x-1),x\in {\Phi }^{+}\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{and}\hfill \\ \hfill & \Re \left(x\right)=\left\{\begin{array}{c}\psi \left(x\right)-\nu \left(x\right)p(x-1),x\in {\Phi }^{-},\hfill \\ \psi \left(x\right),x\in {\Phi }^0,\hfill \\ \psi \left(x\right)-\omega \left(x\right)q(x+1),x\in {\Phi }^{+}\hfill \end{array}\right.\hfill \end{array}$$

with the interval-boundary conditions given in (2).

Letting $\epsilon =0$ yields the corresponding reduced form of (4) as

$${\mathcal{L}}_{\epsilon }{\omega }_0\left(x\right)={\Re }_0\left(x\right),$$

(5)

where

$$\begin{array}{cc}\hfill & {\mathcal{L}}_{\epsilon }{\omega }_0\left(x\right)=\left\{\begin{array}{c}\mu \left(x\right){\omega }_0\left(x\right)+\gamma \left(x\right){\omega }_0(x+1),x\in {\Phi }^{-},\hfill \\ \mu \left(x\right){\omega }_0\left(x\right)+\nu \left(x\right){\omega }_0(x-1)+\gamma \left(x\right){\omega }_0(x+1),x\in {\Phi }^0,\hfill \\ \mu \left(x\right){\omega }_0\left(x\right)+\nu \left(x\right){\omega }_0(x-1),x\in {\Phi }^{+}\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{and}\hfill \\ \hfill & {\Re }_0\left(x\right)=\left\{\begin{array}{c}\psi \left(x\right)-\nu \left(x\right)p_0(x-1),x\in {\Phi }^{-},\hfill \\ \psi \left(x\right),x\in {\Phi }^0,\hfill \\ \psi \left(x\right)-\gamma \left(x\right)q_0(x+1),x\in {\Phi }^{+}\hfill \end{array}\right.\hfill \end{array}$$

In fact the reduced solution ${\omega }_0\left(x\right)$ of (5) need not satisfy ${\omega }_0\left(0\right)=\omega \left(0\right)$ and ${\omega }_0\left(3\right)=\omega \left(3\right)$. As a result, the solution exhibits boundary layers at $x=0$ and $x=3$. Moreover, at the points $x=1$ and $x=2$, the reduced solution may not satisfy $\omega \left(1^{-}\right)=\omega \left(1^{+}\right)$, ${\omega }^{\prime }\left(1^{-}\right)={\omega }^{\prime }\left(1^{+}\right)$, $\omega \left(2^{-}\right)=\omega \left(2^{+}\right)$ and ${\omega }^{\prime }\left(2^{-}\right)={\omega }^{\prime }\left(2^{+}\right)$, where $x^{-}$ and $x^{+}$ denote the left and right sides limiting value of x at 1 and 2, respectively. As a result, the solution involves strong interior layers at $x=1$ and $x=2$ [15].

Lemma 2.1

(Continuous maximum principle). Suppose $\pi \left(x\right)\in \mathtt{M}$ satisfying $\pi \left(x\right)\ge 0$ for $x\in \partial \Phi$ and ${\mathcal{L}}_{\epsilon }\pi \left(x\right)\ge 0$ for $x\in \Phi$. Then $\pi \left(x\right)\ge 0$ for $x\in \overline{\Phi }$.

Proof. 

Let $\lambda \in \overline{\Phi }$ and suppose that $\pi \left(\lambda \right)={min}_{\overline{\Phi }}\pi \left(x\right)<0$. Then, from the given hypothesis, we observe that $\lambda \notin \partial \Phi$, and by the extreme value theorem ${\pi }_x\left(\lambda \right)=0$ and ${\pi }_{xx}\left(\lambda \right)\ge 0$.

Case 1: On the interval $(0,1]$, we have

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }\pi \left(\lambda \right)=& -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+\mu \left(\lambda \right)\pi \left(\lambda \right)+\gamma \left(\lambda \right)\pi (\lambda +1)\hfill \\ \hfill \le & -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+\mu \left(\lambda \right)\pi \left(\lambda \right)+\gamma \left(\lambda \right)\pi \left(\lambda \right),as\gamma \left(\lambda \right)\pi (\lambda +1)\le \gamma \left(\lambda \right)\pi \left(\lambda \right)for\gamma \left(\lambda \right)<0,\hfill \\ \hfill =& -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+[\mu \left(\lambda \right)+\gamma \left(\lambda \right)]\pi \left(\lambda \right)\le 0,as\mu \left(\lambda \right)+\gamma \left(\lambda \right)>0.\hfill \end{array}$$

Case 2: On the interval $(1,2]$, we get

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }\pi \left(\lambda \right)=& -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+\mu \left(\lambda \right)\pi \left(\lambda \right)+\nu \left(\lambda \right)\pi (\lambda -1)+\gamma \left(\lambda \right)\pi (\lambda +1)\le -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+\mu \left(\lambda \right)\pi \left(\lambda \right)\hfill \\ \hfill & +\nu \left(\lambda \right)\pi (\lambda -1)+\gamma \left(\lambda \right)\pi \left(\lambda \right),\hfill \\ \hfill =& -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+[\mu \left(\lambda \right)+\nu \left(\lambda \right)+\gamma \left(\lambda \right)]\pi \left(\lambda \right)+\nu \left(\lambda \right)[\pi (\lambda -1)-\pi \left(\lambda \right)]\hfill \\ \hfill \le & 0,as\mu \left(\lambda \right)+\nu \left(\lambda \right)+\gamma \left(\lambda \right)>0and\pi (\lambda -1)-\pi \left(\lambda \right)\le 0.\hfill \end{array}$$

Case 3: On the interval $(2,3)$, we get

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }\pi \left(\lambda \right)=& -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+\mu \left(\lambda \right)\pi \left(\lambda \right)+\nu \left(\lambda \right)\pi (\lambda -1)\le -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+\mu \left(\lambda \right)\pi \left(\lambda \right)+\nu \left(\lambda \right)\pi (\lambda -1),\hfill \\ \hfill =& -\epsilon \frac{d^2\pi \left(\lambda \right)}{d{x}^2}+[\mu \left(\lambda \right)+\nu \left(\lambda \right)]\pi \left(\lambda \right)+\nu \left(\lambda \right)[\pi (\lambda -1)-\pi \left(\lambda \right)]\hfill \\ \hfill \le & 0,as\mu \left(\lambda \right)+\nu \left(\lambda \right)>0and\pi (\lambda -1)-\pi \left(\lambda \right)\le 0.\hfill \end{array}$$

All the above cases contradict to the hypothesis, which imply that our supposition is wrong. Thus, $\pi \left(x\right)\ge 0$ for $x\in \overline{\Phi }$.    □

Lemma 2.2

(Continuous stability estimate). The continuous solution $\omega \left(x\right)$ of (1)-() is bounded as

$$\left|\omega \left(x\right)\right|\le \frac{\parallel {\mathcal{L}}_{\epsilon }\omega \left(x\right)\parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}$$

Proof. 

Define barrier functions as

$${\varpi }^{\pm }\left(x\right)=\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \left(x\right)\parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\pm \omega \left(x\right).$$

(6)

Then, we have the following results.

$$\begin{array}{cc}\hfill & {\varpi }^{\pm }\left(0\right)=\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \left(x\right)\parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\pm \omega \left(0\right)\ge 0.\hfill \\ \hfill & {\varpi }^{\pm }\left(3\right)=\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \left(x\right)\parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\pm \omega \left(3\right)\ge 0.\hfill \end{array}$$

For $x\in (0,1]$, we have

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }{\varpi }^{\pm }\left(x\right)=& -\epsilon \frac{d^2{\varpi }^{\pm }\left(x\right)}{d{x}^2}+\mu \left(x\right){\varpi }^{\pm }\left(x\right)+\gamma \left(x\right){\varpi }^{\pm }(x+1)\hfill \\ \hfill =& \pm \Re \left(x\right)+[\mu \left(x\right)+\gamma \left(x\right)]\left(\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\right)\hfill \\ \hfill \ge & 2\alpha \left(\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\right)\ge 0.\hfill \end{array}$$

For $x\in (1,2]$, we have

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }{\varpi }^{\pm }\left(x\right)=& -\epsilon \frac{d^2{\varpi }^{\pm }\left(x\right)}{d{x}^2}+\mu \left(x\right){\varpi }^{\pm }\left(x\right)+\nu \left(x\right){\varpi }^{\pm }(x-1)+\gamma \left(x\right){\varpi }^{\pm }(x+1)\hfill \\ \hfill =& \pm \Re \left(x\right)+[\mu \left(x\right)+\nu \left(x\right)+\gamma \left(x\right)]\left(\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\right)\hfill \\ \hfill \ge & 2\alpha \left(\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\right)\ge 0.\hfill \end{array}$$

For $x\in (2,3]$, we have

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }{\varpi }^{\pm }\left(x\right)=& -\epsilon \frac{d^2{\varpi }^{\pm }\left(x\right)}{d{x}^2}+\mu \left(x\right){\varpi }^{\pm }\left(x\right)+\nu \left(x\right){\varpi }^{\pm }(x-1)+\gamma \left(x\right){\varpi }^{\pm }(x+1)\hfill \\ \hfill =& \pm \Re \left(x\right)+[\mu \left(x\right)+\nu \left(x\right)]\left(\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\right)\hfill \\ \hfill \ge & 2\alpha \left(\frac{\parallel {\mathcal{L}}_{\epsilon }\omega \parallel }{\alpha }+\underset{x\in \overline{\Phi }}{max}\left\{\right|\omega \left(0\right)|,|\omega \left(3\right)\left|\right\}\right)\ge 0.\hfill \end{array}$$

Applying Lemma 2.1, we obtain ${\varpi }^{\pm }\left(x\right)\ge 0$, $x\in \overline{\Phi }$, which yields the required stability estimate.    □

Lemma 2.3.

Assuming the conditions in (3) hold true, the ${\ell }^{th}$ order derivative of the solution $\omega \left(x\right),x\in \overline{\Phi }$ of (1)-(1) are bounded as

$$\left|\frac{d^{\ell }\omega \left(x\right)}{d{x}^{\ell }}\right|\le \mathtt{C}(1+{\epsilon }^{-\frac{\ell }{2}}),\ell =0,1,2,3,4.$$

Proof. 

For $\ell =0$, it immediately implies Lemma . For $\ell =1$, consider $x\in (0,1]$ and define a neighborhood ${\aleph }_x=(t,t+\sqrt{\epsilon })$, so that $x\in {\aleph }_x$ and ${\aleph }_x\subset [0,1]$. For some $\varsigma \in {\aleph }_x$, the mean value theorem yields

$$\begin{array}{c}\hfill \left|\frac{d\omega \left(\varsigma \right)}{dx}\right|=\left|\frac{\omega (t+\sqrt{\epsilon })-\omega \left(t\right)}{(t+\sqrt{\epsilon })-t}\right|={\epsilon }^{-\frac{1}{2}}|\omega (t+\sqrt{\epsilon })-\omega \left(t\right)|\le 2{\epsilon }^{-\frac{1}{2}}\parallel \omega \parallel .\end{array}$$

(7)

For $x\in {\aleph }_x$ and $s\in (\varsigma ,x)$, we can have

$$\begin{array}{cc}\hfill \left|\frac{d\omega \left(x\right)}{dx}\right|=& \left|\frac{d\omega \left(\varsigma \right)}{dx}+{\int }_{\varsigma }^x\frac{d^2\omega \left(s\right)}{d{x}^2}ds\right|\hfill \\ \hfill =& \left|\frac{d\omega \left(\varsigma \right)}{dx}+{\epsilon }^{-1}{\int }_{\varsigma }^x[\mu \left(s\right)\omega \left(s\right)+\nu \left(s\right)p(s-1)+\gamma \left(x\right)\omega (s+1)-\psi \left(s\right)]ds\right|\hfill \\ \hfill \le & \left|\frac{d\omega \left(\varsigma \right)}{dx}\right|+{\epsilon }^{-1}[\parallel \mu \omega \parallel +\parallel \nu p\parallel +\parallel \gamma \omega \parallel +\parallel \psi \parallel \left]\right|{\int }_{\varsigma }^{x}ds|\hfill \\ \hfill \le & 2{\epsilon }^{-\frac{1}{2}}\parallel \omega \parallel +{\epsilon }^{-1}[\parallel \mu \omega \parallel +\parallel \nu p\parallel +\parallel \gamma \omega \parallel +\parallel \psi \parallel ]\sqrt{\epsilon },using\left(\right)\hfill \\ \hfill \le & \mathtt{C}{\epsilon }^{-\frac{1}{2}}\le \mathtt{C}(1+{\epsilon }^{-\frac{1}{2}}).\hfill \end{array}$$

For $\ell =2$, rearrangement of (1) becomes

$$\frac{d^2\omega \left(x\right)}{d{x}^2}={\epsilon }^{-1}[\mu \left(x\right)\omega \left(x\right)+\nu \left(x\right)p(x-1)+\gamma \left(x\right)\omega (x+1)-\psi \left(x\right)],$$

(8)

from which we obtain

$$\begin{array}{c}\hfill \left|\frac{d^2\omega \left(x\right)}{d{x}^2}\right|\le {\epsilon }^{-1}[\parallel \mu \omega \parallel +\nu p\parallel +\gamma \omega \parallel +\parallel \psi \parallel ]\le \mathtt{C}{\epsilon }^{-1}\le \mathtt{C}(1+{\epsilon }^{-1}).\end{array}$$

For $\ell =3$, differentiating (8) once becomes

$$\frac{d^3\omega \left(x\right)}{d{x}^3}={\epsilon }^{-1}[{\left(\mu \omega \right)}^{\prime }\left(x\right)+{\left(\nu \left(x\right)p(x-1)\right)}^{\prime }+{\left(\gamma \left(x\right)\omega (x+1)\right)}^{\prime }-{\psi }^{\prime }\left(x\right)],$$

(9)

from which we obtain

$$\left|\frac{d^3\omega \left(x\right)}{d{x}^3}\right|\le {\epsilon }^{-1}[\mathtt{C}+\mathtt{C}{\epsilon }^{-\frac{1}{2}}]\le \mathbb{C}{\epsilon }^{-\frac{3}{2}}\le \mathbb{C}(1+{\epsilon }^{-\frac{3}{2}}).$$

By similar procedure, for $\ell =4$, differentiating (9) yields

$$\left|\frac{d^4\omega \left(x\right)}{d{x}^4}\right|\le {\epsilon }^{-1}[\mathtt{C}+\mathtt{C}{\epsilon }^{-\frac{1}{2}}+\mathtt{C}{\epsilon }^{-1}]\le \mathtt{C}{\epsilon }^{-2}\le \mathtt{C}(1+{\epsilon }^{-2}).$$

Similar procedure hold true on the intervals $(1,2)$ and $(2,3)$, and hence the continuous solution and its derivatives are bounded.    □

3. Formulation of the Numerical Method

Let the interval $[0,3]$ be divided into $3N$ specially constructed sub-intervals, so that the terms containing the shift arguments coincide with the nodal points. Choosing h as an equal mesh size, we have $h=\frac{1}{\xi }=\frac{2}{\eta }=\frac{3}{3N}$, where $\xi$ and $\eta$ are mesh numbers up to the nodal points $x_i=1$ and $x_i=2$, respectively, and satisfying $\xi <\eta <3N$.

Now, rearranging the differential equation in (1) as

$$\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}=\Im \left(x\right),$$

(10)

where $\Im \left(x\right)=\mu \left(x\right)\omega \left(x\right)+\nu \left(x\right)\omega (x-1)+\gamma \left(x\right)\omega (x+1)-\psi \left(x\right)$. At a mesh point $x_i$, we have

$$\epsilon \frac{d^2\omega \left(x_i\right)}{d{x}^2}=\Im \left(x_i\right),$$

(11)

where $\Im \left(x_i\right)=\mu \left(x_i\right)\omega \left(x_i\right)+\nu \left(x_i\right)\omega (x_i-1)+\gamma \left(x_i\right)\omega (x_i+1)-\psi \left(x_i\right)$.

Applying the Nomerov’s finite difference approach [16] on to equation (11), we get

$$\frac{\epsilon }{h^2}({\omega }_{i+1}-2{\omega }_i+{\omega }_{i-1})=\frac{1}{12}({\Im }_{i+1}+10{\Im }_i+{\Im }_{i-1})$$

(12)

To handle the effect of $\epsilon$ in the layers, we introduce the fitting factor $\sigma$ in (12) as

$$\frac{\epsilon \sigma }{h^2}({\omega }_{i+1}-2{\omega }_i+{\omega }_{i-1})=\frac{1}{12}({\Im }_{i+1}+10{\Im }_i+{\Im }_{i-1}),$$

which is expressed as

$$\begin{array}{cc}\hfill \frac{\sigma }{{\rho }^2}({\omega }_{i+1}-2{\omega }_i+{\omega }_{i-1})=& \frac{1}{12}[{\mu }_{i+1}{\omega }_{i+1}+10{\mu }_i{\omega }_i+{\gamma }_{i-1}{\omega }_{i-1}+{\nu }_{i+1}{\omega }_{i+1-N}\hfill \\ \hfill +& 10{\nu }_i{\omega }_{i-N}+{\nu }_{i-1}{\omega }_{i-1-N}]+{\gamma }_{i+1}{\omega }_{i+1+N}+10{\gamma }_i{\omega }_{i+N}+{\gamma }_{i-1}{\omega }_{i-1+N}\hfill \\ \hfill & -\frac{1}{12}[{\psi }_{i+1}+10{\psi }_i+{\psi }_{i-1}],i=0,1,...,3N.\hfill \end{array}$$

(13)

Assuming uniform convergence of the solution and taking $h\to 0$, we obtain

$$\sigma =\frac{{\rho }^2\mu \left(0\right)(2cosh\left(\sqrt{\mu \left(0\right)}\rho \right)+10)}{48{sinh}^2(\sqrt{\mu (j)}\rho /2)},$$

(14)

which is a fitting factor in the left layers for $j=0$ and in the right layers for $j=3$. Using the fitting factor (14), we get

$${\mathcal{L}}_{\epsilon }^{3N}{\omega }_i=\epsilon \sigma D^{+}{D}^{-}{\omega }_i=\frac{1}{12}[{\Im }_{i+1}+10{\Im }_i+{\Im }_{i-1}],$$

which can be re-expressed as

$$\begin{array}{cc}\hfill & \left(\frac{\sigma }{{\rho }^2}-\frac{{\mu }_{i+1}}{12}\right){\omega }_{i+1}-\left(\frac{2\sigma }{{\rho }^2}+\frac{10{\mu }_i}{12}\right){\omega }_i+\left(\frac{\sigma }{{\rho }^2}-\frac{{\mu }_{i-1}}{12}\right){\omega }_{i-1}=\frac{1}{12}[{\nu }_{i+1}{\omega }_{i+1-N}+10{\nu }_i{\omega }_{i-N}\hfill \\ \hfill & +{\nu }_{i-1}{\omega }_{i-1-N}+{\gamma }_{i+1}{\omega }_{i+1+N}+10{\gamma }_i{\omega }_{i+N}+{\gamma }_{i-1}{\omega }_{i-1+N}-({\psi }_{i+1}+10{\psi }_i+{\psi }_{i-1})],\rho =\frac{h}{\sqrt{\epsilon }}.\hfill \end{array}$$

(15)

To accomplish the required discrete scheme, first we set $\epsilon =0$ and determine the reduced problem as

$$\mu \left(x_i\right){\omega }_0\left(x_i\right)+\nu \left(x_i\right){\omega }_0(x_i-1)+\gamma \left(x_i\right){\omega }_0(x_i+1)=\psi \left(x_i\right).$$

(16)

From the reduced scheme (16), compute the approximate values of $\omega (x_i=1)$ and $\omega (x_i=2)$. Then for $i=1,2,\cdots ,N-1$, the scheme in (15) is simplified to

$${\wp }^{+}{\omega }_{i+1}+{\wp }^o{\omega }_i+{\wp }^{-}{\omega }_{i-1}={\Xi }_{i+1}^{-},$$

(17)

with the boundary conditions ${\omega }_0=p\left(0\right)$ and ${\omega }_N={\omega }_0\left(1\right)$, where

$$\begin{array}{cc}\hfill {\wp }^{+}=& \frac{\sigma }{{\rho }^2}-\frac{{\mu }_{i+1}}{12},{\wp }^o=-\frac{2\sigma }{{\rho }^2}-\frac{10{\mu }_i}{12},{\wp }^{-}=\frac{\sigma }{{\rho }^2}-\frac{{\mu }_{i-1}}{12}\hfill \\ \hfill {\Xi }_{i+1}^{-}=& \frac{1}{12}[{\nu }_{i+1}{p}_{i+1-N}+10{\nu }_i{p}_{i-N}+{\nu }_{i-1}{p}_{i-1-N}+{\gamma }_{i+1}{\omega }_{i+1+N}+10{\gamma }_i{\omega }_{i+N}\hfill \\ \hfill & +{\gamma }_{i-1}{\omega }_{i-1+N}-({\psi }_{i+1}+10{\psi }_i+{\psi }_{i-1})],\hfill \end{array}$$

Similarly, for $i=N+1,N+2,\cdots ,2N-1$, the scheme in (15) becomes

$${\wp }^{+}{\omega }_{i+1}+{\wp }^o{\omega }_i+{\wp }^{-}{\omega }_{i-1}={\Xi }_{i+1}^o,$$

(18)

with the boundary conditions ${\omega }_N={\omega }_0\left(1\right)$ and ${\omega }_{2N}={\omega }_0\left(2\right)$ , where

$$\begin{array}{cc}\hfill {\Xi }_{i+1}^o=& \frac{1}{12}[{\nu }_{i+1}{\omega }_{i+1-N}+10{\nu }_i{\omega }_{i-N}+{\nu }_{i-1}{\omega }_{i-1-N}+{\gamma }_{i+1}{\omega }_{i+1+N}+10{\gamma }_i{\omega }_{i+N}+{\gamma }_{i-1}{\omega }_{i-1+N}\hfill \\ \hfill & -({\psi }_{i+1}+10{\psi }_i+{\psi }_{i-1})],\hfill \end{array}$$

and for $i=2N+1,2N+2,...,3N-1$, the scheme in (15) becomes

$${\wp }^{+}{\omega }_{i+1}+{\wp }^o{\omega }_i+{\wp }^{-}{\omega }_{i-1}={\Xi }_{i+1}^{+},$$

(19)

with the boundary conditions ${\omega }_{2N}={\omega }_0\left(2\right)$ and ${\omega }_{3N}=q\left(x_{3N}\right)$ , where

$$\begin{array}{cc}\hfill {\Xi }_{i+1}^{+}=& \frac{1}{12}[{\nu }_{i+1}{\omega }_{i+1-N}+10{\nu }_i{\omega }_{i-N}+{\nu }_{i-1}{\omega }_{i-1-N}+{\gamma }_{i+1}{q}_{i+1+N}+10{\gamma }_i{q}_{i+N}+{\gamma }_{i-1}{q}_{i-1+N}\hfill \\ \hfill & -({\psi }_{i+1}+10{\psi }_i+{\psi }_{i-1})].\hfill \end{array}$$

Thus, we solve the systems obtained in (17)-(19) on a uniform meshes using suitable solver of system of equations following the computational algorithm given below.

Algorithm 1:Computing the numerical solution ${\omega }_{i+1},i=0,1,\cdots ,3N-1$

1.3 1: On $[0,3]$, set a uniform mesh $x_i=ih$, $i=0,1,2,\cdots ,3N$ with $x_0=0$, $x_N=1$, $x_{2N}=2$ and $x_{3N}=3$. 2: Set $\epsilon =0$ and approximate ${\omega }_{i+1}$ from (16) at $x_i=1$ and $x_i=2$. 3: For $i=1,2,\cdots ,3N-1$, 4:       if $i=1,2,\cdots ,N-1$, compute ${\omega }_{i+1}$ from the system (18) 5:       else if $i=N+1,N+2,\cdots ,2N-1$, compute for ${\omega }_{i+1}$ from the system (17) 6:       else compute for ${\omega }_{i+1}$ from the system (19). 7: Simulate the solution ${\omega }_{i+1}$ for all i graphically. 8: Stop

3.1. Uniform Stability and Convergence Analysis

Lemma 3.1

(Discrete maximum principle). Suppose the mesh function ${\vartheta }_i$ be in ${\mathcal{L}}_{\epsilon }^{3N}\left(\overline{\Phi }\right)$ satisfying ${\vartheta }_0\ge 0$, ${\vartheta }_{3N}\ge 0$ and ${\mathcal{L}}_{\epsilon }^{3N}{\vartheta }_i\ge 0$, $i=1,2,\cdots ,3N-1$. Then, ${\vartheta }_i\ge 0$, $i=0,1,\cdots ,3N$.

Proof. 

Let $\xi \in \{0,1,\cdots ,3N\}$ be given such that ${\vartheta }_{\xi }={min}_{0\le i\le 3N}{\vartheta }_i<0$. From the given hypothesis, it is clear that $\xi \notin \{0,3N\}$. Then, we have ${\vartheta }_{\xi }-{\xi }_{\xi -1}\le 0$ and ${\vartheta }_{\xi +1}-{\xi }_{\xi }\ge 0$.

Case 1: For $i=1,\cdots ,N$, we get

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }^{3N}{\vartheta }_{\xi }=& -\epsilon \sigma D^{+}{D}^{-}{\vartheta }_{\xi }+\frac{1}{12}[{\mu }_{\xi +1}{\vartheta }_{\xi +1}+10{\mu }_{\xi }{\vartheta }_{\xi }+{\mu }_{\xi -1}{\vartheta }_{\xi -1}+{\gamma }_{\xi +1}{\vartheta }_{\xi +1+N}\hfill \\ \hfill & +10{\gamma }_{\xi }{\vartheta }_{\xi }+{\gamma }_{\xi -1}{\vartheta }_{\xi -1+N}]\hfill \\ \hfill =& -\frac{\epsilon \sigma }{h^2}[({\vartheta }_{\xi +1}-{\vartheta }_{\xi })-({\vartheta }_{\xi }-{\vartheta }_{\xi -1})]+\frac{1}{12}[{\mu }_{\xi +1}{\vartheta }_{\xi +1}\hfill \\ \hfill & +10{\mu }_{\xi }{\vartheta }_{\xi }+{\mu }_{\xi -1}{\vartheta }_{\xi -1}+{\gamma }_{\xi +1}{\vartheta }_{\xi +1+N}+10{\gamma }_{\xi }{\vartheta }_{\xi }+{\gamma }_{\xi -1}{\vartheta }_{\xi -1+N}]\hfill \\ \hfill \le & -\frac{\epsilon \sigma }{h^2}[({\vartheta }_{\xi +1}-{\vartheta }_{\xi })-({\vartheta }_{\xi }-{\vartheta }_{\xi -1})]+\frac{1}{12}[({\mu }_{\xi +1}+{\gamma }_{\xi +1}){\vartheta }_{\xi +1}+10({\mu }_{\xi }+{\gamma }_{\xi }){\vartheta }_{\xi }\hfill \\ \hfill & +({\mu }_{\xi -1}+{\gamma }_{\xi -1}){\vartheta }_{\xi -1}]\le 0.\hfill \end{array}$$

Case 2: For $i=N+1,...,2N$, we get

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }^{3N}{\vartheta }_{\xi }=& -\epsilon \sigma D^{+}{D}^{-}{\vartheta }_{\xi }+\frac{1}{12}[{\mu }_{\xi +1}{\vartheta }_{\xi +1}+10{\mu }_{\xi }{\vartheta }_{\xi }+{\mu }_{\xi -1}{\vartheta }_{\xi -1}+{\nu }_{\xi +1}{\vartheta }_{\xi +1-N}+10{\nu }_{\xi }{\vartheta }_{\xi -N}\hfill \\ \hfill & +{\nu }_{\xi -1}{\vartheta }_{\xi -1-N}+{\gamma }_{\xi +1}{\vartheta }_{\xi +1+N}+10{\gamma }_{\xi }{\vartheta }_{\xi }+{\gamma }_{\xi -1}{\vartheta }_{\xi -1+N}]\hfill \\ \hfill =& -\frac{\epsilon \sigma }{h^2}[({\vartheta }_{\xi +1}-{\vartheta }_{\xi })-({\vartheta }_{\xi }-{\vartheta }_{\xi -1})]+\frac{1}{12}[{\mu }_{\xi +1}{\vartheta }_{\xi +1}+10{\mu }_{\xi }{\vartheta }_{\xi }+{\mu }_{\xi -1}{\vartheta }_{\xi -1}\hfill \\ \hfill & +{\nu }_{\xi +1}{\vartheta }_{\xi +1-N}+{\nu }_{\xi }{\vartheta }_{\xi -N}+{\nu }_{\xi -1}{\vartheta }_{\xi -1-N}+{\gamma }_{\xi +1}{\vartheta }_{\xi +1+N}+10{\gamma }_{\xi }{\vartheta }_{\xi }+{\gamma }_{\xi -1}{\vartheta }_{\xi -1+N}]\hfill \\ \hfill \le & -\frac{\epsilon \sigma }{h^2}[({\vartheta }_{\xi +1}-{\vartheta }_{\xi })-({\vartheta }_{\xi }-{\vartheta }_{\xi -1})]+\frac{1}{12}[({\mu }_{\xi +1}+{\nu }_{\xi +1}+{\gamma }_{\xi +1}){\vartheta }_{\xi +1}\hfill \\ \hfill & +10({\mu }_{\xi }+{\nu }_{\xi }+{\gamma }_{\xi }){\vartheta }_{\xi }+({\mu }_{\xi -1}+{\nu }_{\xi -1}+{\gamma }_{\xi -1}){\vartheta }_{\xi -1}]\le 0.\hfill \end{array}$$

Case 3: For $i=2N+1,\cdots ,3N-1$, we get

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }^{3N}{\vartheta }_{\xi }=& -\epsilon \sigma D^{+}{D}^{-}{\vartheta }_{\xi }+\frac{1}{12}[{\mu }_{\xi +1}{\vartheta }_{\xi +1}+10{\mu }_{\xi }{\vartheta }_{\xi }+{\mu }_{\xi -1}{\vartheta }_{\xi -1}+{\nu }_{\xi +1}{\vartheta }_{\xi +1-N}\hfill \\ \hfill & +10{\nu }_{\xi }{\vartheta }_{\xi -N}+{\nu }_{\xi -1}{\vartheta }_{\xi -1-N}]\hfill \\ \hfill =& -\frac{\epsilon \sigma }{h^2}[({\vartheta }_{\xi +1}-{\vartheta }_{\xi })-({\vartheta }_{\xi }-{\vartheta }_{\xi -1})]+\frac{1}{12}[{\mu }_{\xi +1}{\vartheta }_{\xi +1}+10{\mu }_{\xi }{\vartheta }_{\xi }+{\mu }_{\xi -1}{\vartheta }_{\xi -1}]\hfill \\ \hfill & +{\nu }_{\xi +1}{\vartheta }_{\xi +1-N}+{\nu }_{\xi }{\vartheta }_{\xi -N}+{\nu }_{\xi -1}{\vartheta }_{\xi -1-N}]\hfill \\ \hfill \le & -\frac{\epsilon \sigma }{h^2}[({\vartheta }_{\xi +1}-{\vartheta }_{\xi })-({\vartheta }_{\xi }-{\vartheta }_{\xi -1})]+\frac{1}{12}[({\mu }_{\xi +1}+{\nu }_{\xi +1}){\vartheta }_{\xi +1}\hfill \\ \hfill & +10({\mu }_{\xi }+{\nu }_{\xi }){\vartheta }_{\xi }+({\mu }_{\xi -1}+{\nu }_{\xi -1}){\vartheta }_{\xi -1}]\le 0.\hfill \end{array}$$

The above three cases contradict to the hypothesis. Thus, the supposition fails, so that ${\vartheta }_{\xi }\ge 0$, $\xi =1,2,\cdots ,3N-1$. Thus, ${\vartheta }_i\ge 0$, $\xi =0,1,\cdots ,3N$.    □

Lemma 3.2.

Under the assumptions in (3), the discrete solution ${\omega }_i,i=0,1,\cdots ,3N$ of (15) is estimated as

$$|{\omega }_i|\le \frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_i\parallel }{\alpha }+\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}$$

Proof. 

Introduce two barrier functions as

$${\varpi }_i^{\pm }=\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_i\parallel }{\alpha }+\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\pm {\omega }_i,i=0,1,\cdots ,3N.$$

Then, we obtain the following results.

For $i=0$, ${\varpi }_0^{\pm }=\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_0\parallel }{\alpha }+{max}_{i\in [0,3N]}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\pm {\omega }_0\ge 0$.

For $i=3N$, ${\varpi }_{3N}^{\pm }=\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{3N}\parallel }{\alpha }+{max}_{i\in [0,3N]}\left\{\right|{\omega }_{3N}|,|{\omega }_{3N}\left|\right\}\pm {\omega }_{3N}\ge 0$.

For $i=1,2,\cdots ,N$, we have

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }^{3N}{\varpi }_i^{\pm }=& -\epsilon \sigma D^{+}{D}^{-}{\varpi }_i^{\pm }+\frac{1}{12}[{\mu }_{i+1}{\varpi }_{i+1}^{\pm }+10{\mu }_i{\varpi }_i^{\pm }+{\mu }_{i-1}{\varpi }_{i-1}^{\pm }+{\gamma }_{i+1}{\varpi }_{i+1+N}^{\pm }\hfill \\ \hfill & +10{\gamma }_i{\varpi }_i^{\pm }+{\gamma }_{i-1}{\varpi }_{i-1+N}^{\pm }]\hfill \\ \hfill =& \pm {\mathcal{L}}_{\epsilon }^{3N}{\omega }_i+\frac{1}{12}\left({\mu }_{i+1}\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{i+1}\parallel }{\alpha }+{\mu }_{i+1}\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left({\mu }_i\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_i\parallel }{\alpha }+{\mu }_i\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left({\mu }_{i-1}\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{i-1}\parallel }{\alpha }+{\mu }_{i-1}\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left({\gamma }_{i+1}\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{i+1}\parallel }{\alpha }+{\gamma }_{i+1}\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left({\gamma }_i\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_i\parallel }{\alpha }+{\gamma }_i\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left({\gamma }_{i-1}\frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{i-1}\parallel }{\alpha }+{\gamma }_{i-1}\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right),\hfill \end{array}$$

which implies that

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\epsilon }^{3N}{\varpi }_i^{\pm }\ge & \frac{1}{12}\left(\alpha \frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{i+1}\parallel }{\alpha }+\alpha \underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left(\alpha \frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_i\parallel }{\alpha }+\alpha \underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right)\hfill \\ \hfill & +\frac{1}{12}\left(\alpha \frac{\parallel {\mathcal{L}}_{\epsilon }^{3N}{\omega }_{i-1}\parallel }{\alpha }+\alpha \underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\right),as{\mu }_i+{\gamma }_i\ge \alpha \hfill \\ \hfill \ge & \frac{\alpha }{12}\underset{i\in [0,3N]}{max}\left\{\right|{\omega }_0|,|{\omega }_{3N}\left|\right\}\ge 0.\hfill \end{array}$$

Applying Lemma 3.1, the result implies that $\varpi \ge 0,i=1,2,\cdots ,N-1$, and consequently it yields the required estimate. Similar procedures hold for $i=N+1,N+2,\cdots ,2N-1$ and $i=2N+1,2N+2,\cdots ,3N-1$.    □

Theorem 3.1.

Let $\omega \left(x_i\right)$ and ${\omega }_i$ be the solution of the prtoblem (1) and the scheme (15), respectively. Then, the truncation error is estimated as

$$|\omega \left(x_i\right)-{\omega }_i|\le \mathtt{C}{N}^{-2},i=0,1,\cdots ,3N,$$

with $\mathtt{C}$ independent of i, h and ε.

Proof. 

From the definitions of truncation error [17], we have

$$\begin{array}{c}\hfill |{\mathcal{L}}_{\epsilon }^{3N}(\omega \left(x_i\right)-{\omega }_i)|=|{\mathcal{L}}_{\epsilon }\omega \left(x_i\right)-{\mathcal{L}}_{\epsilon }^{3N}{\omega }_i|=\left|-\epsilon \frac{d^2\omega \left(x_i\right)}{d{x}^2}+\epsilon \sigma D^{+}{D}^{-}{\omega }_i\right|.\end{array}$$

(20)

Applying Taylor’s series expansion on ${\omega }_{i\pm 1}$, we obtain

$$\begin{array}{cc}\hfill & {\omega }_{i+1}={\omega }_i+h\frac{d\omega \left(x_i\right)}{dx}+\frac{h^2}{2!}\frac{d^2\omega \left(x_i\right)}{d{x}^2}+\frac{h^3}{3!}\frac{d^3\omega \left(x_i\right)}{d{x}^3}+\frac{h^4}{4!}\frac{d^4\omega \left(x_i\right)}{d{x}^4}+O\left(h^5\right),\hfill \\ \hfill & {\omega }_{i-1}={\omega }_i-h\frac{d\omega \left(x_i\right)}{dx}+\frac{h^2}{2!}\frac{d^2\omega \left(x_i\right)}{d{x}^2}-\frac{h^3}{3!}\frac{d^3\omega \left(x_i\right)}{d{x}^3}+\frac{h^4}{4!}\frac{d^4\omega \left(x_i\right)}{d{x}^4}+O\left(h^5\right).\hfill \end{array}$$

Plugging the above expansion into (20) and simplifying gives

$$\begin{array}{cc}\hfill |{\mathcal{L}}_{\epsilon }^{3N}(\omega \left(x_i\right)-{\omega }_i)|\le & \mathtt{C}|\epsilon \sigma -\epsilon |∥\frac{d^2\omega \left(t_i\right)}{d{x}^2}∥+\mathtt{C}\epsilon \sigma h^2∥\frac{d^4\omega \left(t_i\right)}{d{x}^4}∥,t_i\in [x_{i-1},x_{i+1}].\hfill \end{array}$$

As it is estimated in [18], we have $|\epsilon \sigma -\epsilon |\le \mathtt{C}{h}^2$. Applying Taylor’s series expansion on $\sigma$ and using Lemma 2.3, we have

$$|{\mathcal{L}}_{\epsilon }^{3N}(\omega \left(x_i\right)-{\omega }_i)|\le \mathtt{C}{h}^2=\mathtt{C}{N}^{-2},ash=\frac{1}{N}.$$

Thus, by the discrete maximum principle, we get $|\omega \left(x_i\right)-{\omega }_i|\le \mathtt{C}{N}^{-2},i=0,1,\cdots ,3N$.    □

4. Numerical Examples and Discusions

In this section, the developed finite difference scheme is validated by carrying out numerical experiments. The effect of the perturbation parameter and the mixed shift arguments on the solution is determined by computing the point wise maximum absolute errors (MAEs) and convergence rates (CRs). Since computation of the exact solution is difficult, to compute MAEs and CRs we used the double mesh principle [18] as given in Algorithm 2 below.

Algorithm 2:Computations of MAEs and corresponding CRs

1.3 1: From (16)–(19), compute ${\omega }_{i+1}^h$ and ${\omega }_{i+1}^{\frac{h}{2}}$ at each mesh point. 2: For $i=0,1,\cdots ,3N$, compute 3:             $e_{\epsilon }^h=max|{\omega }_{i+1}^h-{\omega }_{i+1}^{\frac{h}{2}}|$, ▹ Point-wise absolute error 4:             $e^h=max\left(e_{\epsilon }^h\right)$, ▹ Uniform absolute error 5:             $r_{\epsilon }^h={log}_2(e_{\epsilon }^h/e_{\epsilon }^{\frac{h}{2}})$, ▹ Point-wise convergence rate 6:             $r^h={log}_2(e^h/e^{\frac{h}{2}})$, ▹ Uniform convergence rate. 7: Sketch log-log plots of MAE versus N.

Example 4.1.

[12] Consider the the SP-FDE given as

$$\begin{array}{cc}\hfill & -\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}+5\omega \left(x\right)-\omega (x-1)-\omega (x+1)=1,x\in \mathtt{M}\hfill \\ \hfill & \left\{\begin{array}{c}\omega \left(x\right)=1,x\in [-1,0],\hfill \\ \omega \left(x\right)=1,x\in [3,4].\hfill \end{array}\right.\hfill \end{array}$$

Example 4.2.

[12] Consider the the SP-FDE given as

$$\begin{array}{cc}\hfill & -\epsilon \frac{d^2\omega \left(x\right)}{d{x}^2}+(5+x^2)\omega \left(x\right)-\omega (x-1)-\omega (x+1)=x,x\in \mathtt{M}\hfill \\ \hfill & \left\{\begin{array}{c}\omega \left(x\right)=1,x\in [-1,0],\hfill \\ \omega \left(x\right)=1,x\in [3,4].\hfill \end{array}\right.\hfill \end{array}$$

The examples are treated applying the formulated numerical method. The error analyses of the solutions are also performed by computing MAEs and CRs by using the double mesh principle as given in Algorithm 2 and the obtained results are displayed in Table 1 and Table 2. From Table 1, one can see that increasing the mesh numbers reduces the error and reducing $\epsilon$ gives a rapidly decreasing MAEs for each mesh number. The maximum errors are shown when $\epsilon =2^{-1}$ with the corresponding second order convergence rates. From Table 2, we observe that if the mesh number is increased, the MAE is decreasing and decreasing $\epsilon$ yields stable maximum absolute errors and corresponding second order convergence rates. The accuracy of the present method is compared with previous published work for both examples as given in Table 3 and Table 4. From these tables, it is shown that the present fitted operator difference scheme is second order convergent method and the fitted mesh method developed in [12] is almost first order convergent method.

The numerical solutions of example 4.1 and 4.2 are shown in Figure 1 and Figure 2, respectively, for different values of $\epsilon$. From the figures, it is shown that decreasing the perturbation parameter reduces the width of the layers with more steep gradients. In Figure 3, log-log plots are given, which show that increasing the number of meshes decreases the MAEs for different value of $\epsilon$.

5. Conclusion and Future Plans

In the present research work, we have developed a parameter-robust numerical for computing SP-FDEs with mixed large shifts. Such equations are a class of differential equations that incorporate shift effects and small perturbation parameters. The presence of the perturbation parameter and large mixed shifts adds complexity to the problem to be solved analytically, and applying standard numerical methods. The proposed numerical scheme tackles the challenges posed by the perturbation parameter and shift arguments yielding accurate approximate solution. We have demonstrated the effectiveness and robustness of the scheme by two numerical experiments, showing that it consistently produces accurate results across a wide range of parameter values. While this research work has provided an efficient and accurate numerical method for the considered singularly perturbed ordinary differential equations with mixed large shifts, there are several avenues for future exploration and enhancement. In the near future, the method is extended to singularly perturbed functional partial differential equations in one and two dimensions.