Shilpkala T. Mane
| Ram Kishun Lodhi*
© 2024 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
In this article, we investigate a hybrid difference scheme for finding the numerical solution of a singularly perturbed second-order reaction-diffusion problem with a discontinuous source term. Such types of problems arise in the modeling of semiconductor devices and geophysical fluid dynamics etc. Solutions of these types of problems are difficult to obtain due to the presence of boundary and interior layers. A hybrid difference scheme i.e., cubic spline method and central finite difference approach, are applied on a fine region and coarse region, respectively. Shishkin mesh is utilized to generate the mesh point for the given domain. We use a second-order hybrid difference operator at the point of discontinuity. The solution rapidly changes in the interior layers and boundary layer. Truncation error is studied, and the stability of the method is analyzed. The proposed method is implemented on two problems, and numerical results are compared with the existing method, which shows that the proposed method is efficient for reducing maximum absolute errors and increasing the rate of convergence.
singularly perturbed problem, reaction-diffusion, discontinuous source term, hybrid difference scheme, cubic spline method, central finite difference approach, Shishkin mesh, interior layer, truncation error
A challenge in the modeling of semiconductor devices [1] is the source of motivation for the knowledge of singularly perturbed differential equations having discontinuous source terms (DST). A differential equation (DE) whose largest derivative term is multiplied by a small positive parameter ε and contains a point of discontinuity in the source term and/or diffusion coefficient is known as a singularly perturbed differential equation with a discontinuous source term. These types of problems have an immense application in applied mathematics and engineering fields such as chemical reactor theory, reaction-diffusion process, quantum mechanics and electrical circuits, etc. These problems frequently have boundary and interior layers where the solution rapidly differs. For the smooth data, the solution of the problem exhibits a steep gradient in some part or all of the boundary of an interval, and for the non-smooth type data, the solution produces extra interior layers. Two main types of singularly perturbed problems are taking attention: problems of the reaction-diffusion kinds where the order is decreased by 1 and convection-diffusion kinds where the order is decreased by 2. In recent years, the singularly perturbed boundary value problem (SPBVP) has drawn much attention from several researchers. Various numerical techniques are implemented in order to obtain the numerical solution of singularly perturbed problems in the literature, more information can be found in the books [2-5]. Here we consider a second-order reaction-diffusion kind boundary value problem:
$L_{\varepsilon} \equiv-\varepsilon u^{\prime \prime}(t)+b(t) u(t)=g(t), \quad t \in \Psi^{-} \cup \Psi^{+}$ (1)
$u(0)=u_0, u(1)=u_1$ (2)
where, $0<\varepsilon \ll 1$ is a singular perturbation parameter, $b(t) \geq \beta(t)>0$ and $g(t)$ are sufficiently smooth on $\bar{\Psi} \backslash\{w\}$. Here single discontinuity in the source term has occurred at the point $w \in \Psi$, where $\Psi=(0,1), \Psi^{-}=(0, w), \Psi^{+}=(w, 1)$. Since g is discontinuous at w, the solution u(t) does not have a continuous second-order derivative at w , i.e., $u \notin C^2(\Psi)$. We denote the jump at w in any function with $[g](w)=g(w+)-g(w-)$.
Farrell et al. [6] applied a finite difference method (FDM) using Shishkin mesh for a second-order SPBVP having non-smooth data and obtained a first-order uniform convergence. Roos and Zarin [7] have solved the same problem using Galerkin finite element method using the Backhvalov-Shishkin mesh with convergence two. Shanthi et al. [8] obtained a numerical solution of second-order SPBVP with discontinuity term in the right-hand side of the inhomogeneous equation containing two parameters using the fitted mesh method and obtained the error bounds. Researchers have solved the SPBVP with discontinuous source terms by numerical techniques [9-12] to achieve a better numerical solution. Cen et al. [13] obtained fourth-order convergence using a high-order finite difference method having interior layers. The numerical solution of the two-point boundary value problem having smooth data using splines has been studied by several researchers [14-22]. Splines are easy to apply with many benefits, as once the solution is calculated, one can easily get the spline interpolation information at any nodal point. For more details about the splines, reader may refer to the books in references [23-25]. Motivated by the above work of different researchers, we developed a hybrid method of second-order SPBVP having DST and discontinuous reaction coefficients. The method consists of a standard central FDM applied in the outer layer, and the cubic spline method (CSM) is used in the inner layer to preserve the monotonicity property. This method minimizes the absolute error as compared to the other existing method in the literature for the considered problem.
The paper is organized as follows: Section 2, contains the maximum principle, existence theorem, stability, and bounds on solution and derivatives. In Section 3, we describe the derivation of the cubic spline method, hybrid scheme, and Shishkin mesh. Error analysis has been addressed in Section 4. Numerical results are provided in Section 5, and the final conclusion is summarized in Section 6.
This section provides the derivation of the existence theorem of a given problem. Also, derivation of maximum principle is provided, which is important and useful for derivation of stability result.
Theorem 1: The SPBVP (1-2) has a solution $u \in C^1(\Psi) \cap C^2\left(\Psi^{-} \cup \Psi^{+}\right)$.
Proof: Proof by construction. Let $S_1, S_2$ be the particular solution of DE: $\begin{aligned} & -\varepsilon S_1+b(t) S_1=g, t \in \Psi^{-} \text {and }-\varepsilon S_2+b(t) S_2=g, t \in \Psi^{+} \text {. }\end{aligned}$.
Consider the function:
$S(t)=\left(\begin{array}{l}S_1(t)+\left(u(0)-S_1(0)\right) \xi_1(t)+E \xi_2(t), t \in \Psi^{-} \\ S_2(t)+\left(u(1)-S_2(1)\right) \xi_2(t)+F \xi_1(t), t \in \Psi^{+}\end{array}\right.$,
where, $\xi_1(t), \xi_2(t)$ are solutions of BVP given as:
$\begin{array}{ll}-\varepsilon \xi_1^{\prime \prime}+b(t) \xi_1=0, & t \in \Psi, \xi_1(0)=1, \xi_1(1)=0 \\ -\varepsilon \xi_2^{\prime \prime}+b(t) \xi_2=0, & t \in \Psi, \xi_2(0)=1, \xi_2(1)=1\end{array}$
and E, F are constant to be chosen so that $S \in C^1(\Psi)$. Also, note that on $(0,1), 0<\xi_j<1, j=1,2$.
Thus $\xi_1, \xi_2$ cannot have an interval maximum or minimum, and hence $\xi_1^{\prime}<0, \xi_2^{\prime}>0, t \in(0,1)$.
We wish to choose the constant E, F so that $S \in C^1(\Psi)$.
So, we impose the conditions $S\left(w^{-}\right)=S\left(w^{+}\right)$ and $S^{\prime}\left(w^{-}\right)=S^{\prime}\left(w^{+}\right)$.
For the constant E, F to exist, we require:
$\left|\begin{array}{l}\xi_2(w)-\xi_1(w) \\ \xi_2^{\prime}(w)-\xi_1^{\prime}(w)\end{array}\right| \neq 0$.
This follows that $\xi_2^{\prime}(w) \xi_1(w)-\xi_2(w) \xi_1^{\prime}(w)>0$.
Let $L_{\varepsilon}$ be the differential operator, which is defined as:
$L_{\varepsilon} u=-\varepsilon u^{\prime \prime}+b u$.
Then $L_{\varepsilon}$ satisfies the following maximum principle on $\bar{\Psi}$.
Theorem 2: (Maximum Principle) Let a function $u(t) \in C^0(\bar{\Psi}) \cap C^1(\Psi) \cap C^2\left(\Psi^{-} \cup \Psi^{+}\right)$ satisfies, $\begin{array}{r}u(0) \geq 0, u(1) \geq 0, L_{\varepsilon} u(t) \geq 0, \forall t \in \Psi^{-} \cup \Psi^{+},[u] w =0,\left[u^{\prime}\right] w \leq 0 . \text { Then, } u(t) \geq 0, \forall t \in \Psi\end{array}$.
Proof: For the proof refer [9].
An immediate consequence of the maximum principle is the following stability result.
Theorem 3: (Stability Result) Let $u(t)$ be a solution ($L_{\varepsilon}$) then, $\|u(t)\|_{\bar{\Psi}} \leq \max \left\{|u(0)|,|u(1)|, \frac{1}{\beta}\|g\|_{\Psi^{-} u \Psi^{+}}\right\}$.
Proof: Let $\Phi_{+}(t)=R \pm u(t)$, where $R=\max \left\{|u(0)|,|u(1)|, \frac{1}{\beta}\|g\|_{\Psi^{-} u \Psi^{+}}\right\}$.
Obviously, $\Phi_{+}(0) \geq 0$, $\Phi_{+}(1) \geq 0$ and for each $t \in \Psi^{-} \cup \Psi^{+}$, $L_{\varepsilon} \Phi_{+}(t)=b(t) R \pm L_{\varepsilon} u(t) \geq \beta R \pm g(t) \geq 0$.
Furthermore, since $u \in C^{\prime}(\Psi)$, $\Phi_{+}(w)= \pm[u](w)=0$ and $\Phi_{+}^{\prime}(w)= \pm\left[u^{\prime}\right](w)=0$.
It follows from the maximum principle that: $\Phi_{+}(t) \geq 0 \quad \forall t \in \bar{\Psi}$, which gives the desired bound on $u$.
To derive the parameter-robust properties of the numerical method, we decompose the solution $u_{\varepsilon}$ of Eqs. (1)-(2) into two components i.e. regular components $u_{1 \varepsilon}$ and singular components $u_{2 \varepsilon}$. We write solution $u_{\varepsilon}$ as $u_{\varepsilon}=u_{1 \varepsilon}(t)+u_{2 \varepsilon}(t)$.
The regular component $u_{1 \varepsilon}$ is defined to be solution of an inhomogeneous system: $L_{\varepsilon} u_{1 \varepsilon}=g$ on $\Psi^{-} \cup \Psi^{+}$ with conditions $\begin{gathered}u_{1 \varepsilon}(0)=\frac{g(0)}{b(0)}, u_{1 \varepsilon}(w-)=\frac{g(w-)}{b(w)}, u_{1 \varepsilon}(w+)=\frac{g(w+)}{b(w)}, u_{1 \varepsilon}(1)=\frac{g(1)}{b(1)}\end{gathered}$.
The singular component $u_{2 \varepsilon}$ is defined to be solution of homogeneous system $L_{\varepsilon} u_{2 \varepsilon}=0$ on $\Psi^{-} \cup \Psi^{+}$ which satifies $\begin{aligned} & {\left[u_{2 \varepsilon}(w)\right]=-\left[u_{1 \varepsilon}(w)\right],\left[u_{2 \varepsilon}(w)\right]=-\left[u_{1 \varepsilon}(w)\right], u_{2 \varepsilon}(0)=} \\ & u_{\varepsilon}(0)-u_{1 \varepsilon}(0), u_{2 \varepsilon}(1)=u_{\varepsilon}(1)-u_{1 \varepsilon}(1),\end{aligned}$ where $[u](w)=u(w+)-u(w-)$ denotes the jump of any function at discontinuity point w.
Theorem 4: For each integer i, satisfies 0≤i≤4, the smooth $u_{1 \varepsilon}(t)$ and singular $u_{2 \varepsilon}(t)$ satisfy the bounds.
$\left|u_{1 \varepsilon}^i(t)\right| \leq\left\{\begin{array}{l}C\left(1+\varepsilon^{1-i / 2} \lambda_1(t)\right), t \in \Psi^{-} \\ C\left(1+\varepsilon^{1-i / 2} \lambda_2(t)\right), t \in \Psi^{+}\end{array}\right.$
$\left|u_{2 \varepsilon}^i(t)\right| \leq\left\{\begin{array}{l}C\left(\varepsilon^{-i / 2} \lambda_1(t)\right), t \in \Psi^{-} \\ C\left(\varepsilon^{-i / 2} \lambda_2(t)\right), t \in \Psi^{+}\end{array}\right.$
where, C is a constant independent of $\varepsilon$ and $\lambda_1(t)=e^{-t \sqrt{\frac{\beta}{\varepsilon}}}+e^{-(d-t) \sqrt{\frac{\beta}{\varepsilon}}}, \lambda_2(t)=e^{-(t-d) \sqrt{\frac{\beta}{\varepsilon}}}+e^{-(1-t) \sqrt{\frac{\beta}{\varepsilon}}}$.
Proof: The bounds on the regular and singular components and their derivatives can be obtained by adopting the procedure fused in Tamilselvan et al. [26].
In this section, we describe CSM, hybrid scheme and Shishkin mesh to generate mesh points.
3.1 Cubic spline difference scheme
Cubic spline is a continuous function and has continuous first and second derivatives. The basic interpolation and approximate integration procedures both benefit from the use of the spline as a tool. Its efficiency in numerical differentiation, though, stands out as a standout feature. For SPBVP, cubic splines are widely utilized to determine the solution. We derive the cubic spline scheme on the interval $\bar{\Psi}=[0,1]$, where the step size is defined by $h_i=t_{i+1}-t_i, i=0, \ldots, N-1$. Following is the main properties of cubic spline interpolating polynomials $S(t)$:
(i) $S(t)$ is a polynomial of degree three on each subinterval $\left[t_i, t_{i+1}\right], 0 \leq i \leq N-1$.
(ii) $S(t)$ is twice continuously differentiable function on $0 \leq t \leq 1$.
(iii) $S\left(t_i\right)=u\left(t_i\right), 0 \leq i \leq N-1$.
The expression of cubic spline is as follows:
$\begin{aligned} & S(t)=\frac{\left(t_{i+1}-t\right)^3}{6 h_i} M_i+\frac{\left(t-t_i\right)^3}{6 h_i} M_{i+1}+\left(u\left(t_i\right)-\frac{h_i^2}{6} M_i\right)\left(\frac{t_{i+1}-t}{h_i}\right) +\left(u\left(t_{i+1}\right)-\frac{h_i^2}{6} M_{i+1}\right)\left(\frac{t-t_i}{h_i}\right), \\ & \mathrm{t}_i \leq t \leq t_{i+1}, \\ & i=0, \mathrm{~K}, N-1 .\end{aligned}$ (3)
where, $M_i=S^{\prime \prime}\left(t_i\right), i=0, \ldots, N$.
The first derivative of $S(t)$ is given by:
$\begin{aligned} & S^{\prime}(t)=-\frac{\left(t_{i+1}-t\right)^2}{2 h_i} M_i+\frac{\left(t-t_i\right)^2}{2 h_i} M_{i+1} +\left(\frac{u\left(t_{i+1}\right)-u\left(t_i\right)}{h_i}\right)-\left(\frac{M_{i+1}-M_i}{6}\right) h_i, \\ & \mathrm{t}_i \leq t \leq t_{i+1}, \quad i=0, \mathrm{~K}, N-1 .\end{aligned}$ (4)
and the second derivative of $S(t)$ is:
$S^{\prime \prime}(t)=\frac{\left(t_{i+1}-t\right)}{h_i} M_i+\frac{\left(t-t_i\right)}{h_i} M_{i+1}$. (5)
From Eq. (4), left-hand first derivative and right-hand first derivative is given by:
$\begin{aligned} & S^{\prime}\left(t_i-\right)=\frac{h_{i-1}}{6} M_{i-1}+\frac{h_{i-1}}{3} \\ & M_i+\frac{u\left(t_i\right)-u\left(t_{i-1}\right)}{h_i} .\end{aligned}$ (6)
and
$\begin{aligned} & S^{\prime}\left(t_i+\right)=-\frac{h_i}{3} M_i-\frac{h_i}{6} M_{i+1} +\frac{u\left(t_{i+1}\right)-u\left(t_i\right)}{h_i} .\end{aligned}$ (7)
From Eq. (3) and Eq. (5), the function $S(t)$ and $S^{\prime \prime}(t)$ are continuous on $\bar{\Psi}$ and for $S^{\prime}(t)$ to be continuous at the interior points $t_i$, we have from Eqs. (6)-(7), the following famous continuity condition:
$\begin{aligned} & \frac{h_{i-1}}{6} M_{i-1}+\left(\frac{h_i+h_{i-1}}{3}\right) M_i+\frac{h_i}{6} M_{i+1} =\left(\frac{u\left(t_{i+1}\right)-u\left(t_i\right)}{h_i}\right)-\left(\frac{u\left(t_i\right)-u\left(t_{i-1}\right)}{h_{i-1}}\right), \\ & i=1, \mathrm{~K}, N-1 .\end{aligned}$ (8)
Substituting $-\varepsilon M_j+c\left(t_j\right) u\left(t_j\right)=f\left(t_j\right), j=i, i \pm 1$, in Eq. (8), following system is obtained:
$\begin{aligned} & {\left[\frac{-3 \varepsilon}{h_{i-1}\left(h_i+h_{i-1}\right)}+\frac{h_{i-1}}{2\left(h_i+h_{i-1}\right)} c_{i-1}\right] u_{i-1}+}{\left[\frac{3 \varepsilon}{h_{i-1} h_i}+c_i\right] u_i+\left[\frac{-3 \varepsilon}{h_i\left(h_i+h_{i-1}\right)}+\frac{h_i}{2\left(h_i+h_{i-1}\right)} c_{i+1}\right] u_{i+1}} \\ & =\left[\frac{h_{i-1}}{2\left(h_i+h_{i-1}\right)}\right] f_{i-1}+f_i+\left[\frac{h_i}{2\left(h_i+h_{i-1}\right)}\right] f_{i+1} .\end{aligned}$ (9)
Using the boundary conditions we obtain the system of linear algebraic equations, which gives the approximate solution i.e. $u_0, u_1, \ldots, u_N$ at nodal points $t_0, t_1, \ldots, t_N$.
3.2 Hybrid scheme
The problem (1-2) is discretized by using the CSM and central FDM on fine region and outer regions respectively. A second order one sided difference approximation is used at the point of discontinuity $t_i=t_{\frac{N}{2}}=w$, given by:
$\begin{gathered}L_t^N u_{\frac{N}{2}}=\frac{-u_{\frac{N}{2}+2}+4 u_{\frac{N}{2}+1}-3 u_{\frac{N}{2}}}{2 h} \\ -\frac{u_{\frac{N}{2}-2}-4 u_{\frac{N}{2}-1}+3 u_{\frac{N}{2}}}{2 h}=0\end{gathered}$ (10)
Hence hybrid scheme is given by:
$L^N u_i=\left\{\begin{array}{lr}L_p^N u_i, & i \neq \frac{N}{2} \\ L_i^N u_i, & i=\frac{N}{2}\end{array}\right.$ (11)
$u(0)=u_0, u(1)=u_1$ (12)
$\begin{aligned} & L_p^N u_i \equiv t_i^{-} u_{i-1}+t_i^c u_i+t_i^{+} u_{i+1} =s_i^{-} f_{i-1}+s_i^c f_i+s_i^{+} f_{i+1}, \\ & i=1, \mathrm{~K}, \frac{N}{2}-1, \frac{N}{2}+1, \mathrm{~K}, N-1 .\end{aligned}$ (13)
where,
$\left\{\begin{array}{l}t_i^{-}=\frac{-3 \varepsilon}{h_{i-1}\left(h_i+h_{i-1}\right)}+\frac{h_{i-1}}{\left(h_i+h_{i-1}\right)} c_{i-1} ; \\ t_i^c=\frac{3 \varepsilon}{h_{i-1} h_i}+c_i ; \\ t_i^{+}=\frac{-3 \varepsilon}{h_i\left(h_i+h_{i-1}\right)}+\frac{h_i}{2\left(h_i+h_{i-1}\right)} c_{i+1} ; \\ s_i^{-}=\frac{h_{i-1}}{2\left(h_i+h_{i-1}\right)} ; s_i^c=1 ; \\ s_i^{+}=\frac{h_i}{2\left(h_i+h_{i-1}\right)} ; \\ \text { for } i=1, \mathrm{~K}, \frac{N}{8}-1, \frac{3 N}{8}+1, \mathrm{~K}, \frac{N}{2}-1, \\ \frac{N}{2}+1, \mathrm{~K}, \frac{5 N}{8}-1, \frac{7 N}{8}+1, \mathrm{~K}, N-1 .\end{array}\right.$ (14)
and
$\left\{\begin{array}{l}t_i^{-}=\frac{-2 \varepsilon}{h_{i-1}\left(h_i+h_{i-1}\right)} ; \\ t_i^c=\frac{2 \varepsilon}{h_{i-1}\left(h_i+h_{i-1}\right)}+\frac{2 \varepsilon}{h_i\left(h_i+h_{i-1}\right)}+c_i \\ t_i^{+}=\frac{-2 \varepsilon}{h_i\left(h_i+h_{i-1}\right)} ; \\ s_i^{-}=0 ; s_i^c=1 ; s_i^{+}=0 ; \\ \text { for } i=\frac{N}{8}, \mathrm{~K}, \frac{3 N}{8}, \frac{5 N}{8}, \mathrm{~K}, \frac{7 N}{8} .\end{array}\right.$ (15)
with given boundary conditions.
The associated matrix for Eqs. (11)-(13) is not an M-matrix. To hold the monotonicity property of matrix for the hybrid scheme, we first find the value of $u_{\frac{N}{2}-2}, u_{\frac{N}{2}+2}$ from Eq. (13) by using CSM.
$u_{\frac{N}{2}-2}=\frac{1}{t_{\frac{N}{2}-1}^{-}}\left[\begin{array}{l}-t_{\frac{N}{2}-1}^c u_{\frac{N}{2}-1}-t_{\frac{N}{2}-1}^{+} u_{\frac{N}{2}}+s_{\frac{N}{2}-1}^{-} f_{\frac{N}{2}-2} \\ +s_{\frac{N}{2}-1}^c f_{\frac{N}{2}-1}+s_{\frac{N}{2}-1}^{+} f_{\frac{N}{2}}\end{array}\right]$
and
$u_{\frac{N}{2}+2}=\frac{1}{t_{\frac{N}{2}+1}^{+}}\left[\begin{array}{c}-t_{\frac{N}{2}+1}^{-} u_{\frac{N}{2}}-t_{\frac{N}{2}+1}^c u_{\frac{N}{2}+1}+s_{\frac{N}{2}+1}^{-} f_{\frac{N}{2}} \\ +s_{\frac{N}{2}+1}^c f_{\frac{N}{2}+1}+s_{\frac{N}{2}+1}^{+} f_{\frac{N}{2}+2}\end{array}\right]$
Substituting the values of $u_{\frac{N}{2}-2}, u_{\frac{N}{2}+2}$ in $L_t^N$. We get,
$\begin{aligned} & L_S^N u_i=\left(4+\frac{t_N^c}{t_{\frac{N}{2}-1}^{-}}\right) u_{\frac{N}{2}-1}+\left(-6+\frac{t_{\frac{N}{2}+1}^{-}}{t_{\frac{N}{2}+1}^{+}}+\frac{t_{\frac{N}{2}-1}^{-}}{t_{\frac{N}{2}-1}^{+}}\right) u_{\frac{N}{2}}+ \left(4+\frac{\frac{t_N^c+1}{2}}{t_{\frac{N}{2}+1}^{-}}\right) u_{\frac{N}{2}+1}=\frac{s_{\frac{N}{2}+1}^{-}}{t_{\frac{N}{2}+1}^{+}} f_{\frac{N}{2}}+\frac{\frac{s_N^c+1}{t^{+}}}{\frac{t_{+1}}{2}} f_{\frac{N}{2}+1}+\frac{s_{\frac{N}{2}+1}^{+}}{t_{\frac{N}{2}+1}^{+}} f_{\frac{N}{2}+2} \end{aligned}$$+\frac{\frac{s_N^{-}}{2}-1}{t_{\frac{N}{2}-1}^{-}} f_{\frac{N}{2}-2}+\frac{s_{\frac{N}{2}-1}^c}{t_{\frac{N}{2}-1}^{-}} f_{\frac{N}{2}-1}+\frac{\frac{s_N^{+}}{2}-1}{t_{\frac{N}{2}-1}^{-}} f_{\frac{N}{2}}$
A hybrid difference approach that is more precise is as follows:
$L^N u_i= \begin{cases}L_p^N u_i, & i \neq \frac{N}{2} \\ L_S^N u_i, & i=\frac{N}{2}\end{cases}$ (16)
$u(0)=u_0, \quad u(1)=u_1$. (17)
3.3 Piecewise uniform Shishkin mesh
We use fitted mesh method to construct a mesh on $\Psi^{-} \cup \Psi^{+}$. We divide the intervals $\Psi^{-}$and $\Psi^{+}$ into three sub-intervals:
$\left.\left.0, \mu_1\right), \mu_1, w-\mu_1\right)$ and $\left.w-\mu_1, w\right)$ are sub-interval of $\Psi^{-}$ and $\mu_1$ satisfies $0<\mu_1 \leq \frac{w}{4}$. Now, $w, w+\mu_2, w+\mu_2, 1-\mu_2$ and $1-\mu_2$ are sub-interval of $\Psi^{+}$ and µ2 satisfies $0<\mu_2 \leq \frac{1-w}{4}$.
We generate $\frac{N}{8}$ mesh points on $\left.\left.0, \mu_1\right), w-\mu_1, w\right), w, w+\mu_2, 1-\mu_2, 1$ and $\frac{N}{4}$ mesh points on $\left.\mu_1, w-\mu_1\right), w+\mu_2, 1-\mu_2$. The interior points of the mesh are given as:
$\Psi_{\varepsilon}^N=\left\{t_i, 1 \leq i \leq \frac{N}{2}-1\right\} \cup\left\{t_i, \frac{N}{2}+1 \leq i \leq N-1\right\}$ also $t_{\frac{N}{2}}=w$ and $\bar{\Psi}_{\varepsilon}^N=\left\{t_i\right\}_0^N$. For the value of $\mu_1=\frac{w}{4}$ and $\mu_2=\frac{1-w}{4}$ the mesh becomes uniform. The choice of the $\mu_1$ and $\mu_2$ is very crucial for the convergence of the method and this transition points depends on the value of N., $\varepsilon$ which is defined as follows:
$\mu_1=\min \left\{\frac{w}{4}, 2 \sqrt{\frac{\varepsilon}{\beta}} \ln N\right\}$ and $\mu_2=\min \left\{\frac{1-w}{4}, 2 \sqrt{\frac{\varepsilon}{\beta}} \ln N\right\}$
The step size in the intervals $\left.\left.0, \mu_1\right), w-\mu_1, w\right)$ is $h_1=\frac{8 \mu_1}{N}$, in the interval $\left.\mu_1, w-\mu_1\right)$ step size is $h_2=\frac{4\left(w-2 \mu_1\right)}{N}$, in the interval $w, w+\mu_2$ and $1-\mu_2$, 1 step size is $h_3=\frac{8 \mu_2}{N}$ and in the interval $w+\mu_2, 1-\mu_2$ step size is $h_4=\frac{4\left(1-w-2 \mu_2\right)}{N}$.
This section elaborates truncation error (TE) and stability of the proposed method. TE [14] for the hybrid difference scheme is defined as follows:
$T_{i, u}=\frac{\varepsilon h^2}{8} u^{i v}\left(t_i\right)+O\left(N^{-3}\right)$ (18)
for $\begin{aligned} & i=1, \ldots, \frac{N}{8}-1, \frac{3 N}{8}+1, \ldots, \frac{N}{2}-1, \frac{N}{2}+1, \ldots, \frac{5 N}{8}-1, \frac{7 N}{8}+1 \ldots, N .\end{aligned}$.
TE is given as:
$\begin{aligned} & T_{i, u}=\frac{-\varepsilon\left(h_i-h_{i-1}\right)}{3} u^{\prime \prime \prime}\left(t_i\right)+\frac{\varepsilon}{12}\left(\frac{h_i^3+h_{i-1}^3}{h_i+h_{i-1}}\right) u^{i v}\left(t_i\right)+O\left(N^{-3}\right)\end{aligned}$ (19)
In outer region the truncation error for points other than transition points can also be presented as:
$T_{i, u}=\frac{2 \varepsilon}{h_i+h_{i-1}}\left(\frac{Q_3\left(t_i, t_{i-1}, u\right)}{h_{i-1}}+\frac{Q_3\left(t_i, t_{i+1}, u\right)}{h_i}\right)$ (20)
where, $Q_n\left(q_1, q_2, r\right)=\frac{1}{n} \int_{q_1}^{q_2}\left(q_2-\zeta\right) r^{n+1}(\zeta) d \zeta$ is remainder term at the transition points $i=\frac{N}{8}, \frac{3 N}{8}, \frac{5 N}{8}, \frac{7 N}{8}$. The truncation error is given as:
$T_{i, u}=\frac{2 \varepsilon}{h_i+h_{i-1}}\left(\frac{Q_2\left(t_i, t_{i-1}, u\right)}{h_{i-1}}+\frac{Q_2\left(t_i, t_{i+1}, u\right)}{h_i}\right)$ (21)
Theorem 5: Let $u(t)$ and $u_i(t)$ be the solution of Eqs. (1)-(2) and Eqs. (11)-(12), respectively. Then the truncation error satisfies the following bounds:
(i) $\begin{aligned} & \left|T_{i, u}\right| \leq C N^{-1} \varepsilon^{\frac{-1}{2}} \beta^{\frac{-1}{2}} \ln N \text { for } 0 \leq i<\frac{N}{8}, \frac{3 N}{8}<i<\frac{N}{2}, \frac{N}{2}+ 1 & <i<\frac{5 N}{8}, \frac{7 N}{8}<i \leq N .\end{aligned}$
(ii) $\begin{aligned} & \quad\left|T_{i, u}\right| \leq C N^{-1}\left(1-\beta^{\frac{-1}{2}} \varepsilon^{\frac{1}{2}} \ln N\right) \text { for } \frac{N}{8}<i<\frac{3 N}{8}, \frac{5 N}{8}<i< \frac{7 N}{8} .\end{aligned}$
(iii) $\left|T_{i, u}\right| \leq C \varepsilon$ for $i=\frac{N}{8}, \frac{3 N}{8}, \frac{5 N}{8}, \frac{7 N}{8}$.
(iv) (iv) $\left|T_{i, u}\right| \leq C N^{-1} \beta^{\frac{-1}{2}} \ln N$ for $i=\frac{N}{2}$.
Proof: We divide the following cases:
Case (i) For $\left.t_i \in 0, \mu_1\right),\left(w-\mu_1, w\right),\left(w, w+\mu_2\right), 1-\mu_2, 1$.
$h_1=h_3=h^{\prime}$.
From Eq. (18) and Theorem 4 i.e., bounds on the derivatives, we get:
$\left|T_{i, u}\right| \leq \frac{2 \varepsilon^{-1} h C_1\left(\lambda_1+\lambda_2\right)}{8}$.
Using $h^{\prime}=16 N^{-1} \varepsilon^{\frac{1}{2}} \beta^{\frac{-1}{2}} \ln N$ and bounding the exponential functions by constants, it is easy to show that, $\begin{aligned} & \left|T_{i, u}\right| \leq C N^{-1} \varepsilon^{\frac{-1}{2}} \beta^{\frac{-1}{2}} \ln N \text { for } 0 \leq i<\frac{N}{8}, \frac{3 N}{8}<i<\frac{N}{2}, \quad \frac{N}{2}+ 1<i<\frac{5 N}{8}, \frac{7 N}{8}<i \leq N .\end{aligned}$
Case (ii) For $t_i \in\left(\mu_1, w-\mu_1\right)$
Using $h_2=\frac{4\left(w-2 \mu_1\right)}{N}$ and from Eq. (19) and Theorem 4 we get,
$\begin{aligned} & \left|T_{i, u}\right| \leq \frac{C_1 h_2\left(\varepsilon+\mu_1\right)}{12} \leq \frac{C_1 h_2}{12} \leq C N^{-1}\left(w-4 \varepsilon^{1 / 2} \beta^{-1 / 2} \ln N\right) \leq C N^{-1}\left(1-\varepsilon^{1 / 2} \beta^{-1 / 2} \ln N\right) \text { for } \frac{N}{8}<i<\frac{3 N}{8} .\end{aligned}$
Similarly, we can show the bound-on truncation error for $\frac{5 N}{8}<i<\frac{7 N}{8}$.
Case (iii) Transition point at $t_{\frac{N}{8}}=\mu_1$ and from the truncation error in Eq. (21), we have:
$\begin{aligned} & \left|\frac{T}{8}, u\right|=\left|\frac{2 \varepsilon}{h_{i-1}+h_i}\left(\frac{Q_2\left(t_{i-1,} t_i, u\right)}{h_{i-1}}+\frac{Q_2\left(t_{i,} t_{i+1}, u\right)}{h_i}\right)\right| \\ & \leq\left|\frac{2 \varepsilon}{h_1+h_2}\left(\frac{1}{2 h_1} \int_{\frac{N}{8}-1}^{\frac{t_N}{8}}\left(t_{\frac{N}{8}}-\zeta\right) u^{\prime \prime \prime}\left(t_{\frac{N}{8}}\right) d \zeta\right)\right| +\left|\frac{2 \varepsilon}{h_1+h_2}\left(\frac{1}{2 h_2} \int_{\frac{t^{\frac{N}{8}}}{8}}^{\frac{t_N}{8}}\left(t_{\frac{N}{8}+1}-\zeta\right) u^{\prime \prime}\left(t_{\frac{N}{8}+1}\right) d \zeta\right)\right| \\ & \end{aligned}$
Using the bounds on the derivatives and simplifying, we get, $\left|T_{\frac{N}{8}, u}\right| \leq C \varepsilon$.
In the similar way at the remaining transition points, we can obtain the required truncation error.
Case (iv) At the point of discontinuity i.e., $\frac{t_N}{2}=w$.
Using the Eq. (10) at the point of discontinuity and using Taylor’s Theorem, we have:
$\begin{aligned} & \left.\left|T_{\frac{N}{2}, u}\right|=\left\lvert\,-\left[t_{\frac{N}{2}}\right)+\left(h_{\frac{N}{2}}+h_{\frac{N}{2}+1}\right) u^{\prime}\left(t_{\frac{N}{2}}\right)+\left(\frac{h_{N / 2}^2}{2}+h_{\frac{N}{2}} h_{\frac{N}{2}+1}\right) u^{\prime \prime}\left(t_{\frac{N}{2}}\right)+\mathrm{K}\right.\right]+ \\ & 4\left[u\left(t_{\frac{N}{2}}\right)+h_{\frac{N}{2}} u^{\prime}\left(t_{\frac{N}{2}}\right)+\frac{h_N^2 / 2}{2} u^{\prime \prime}\left(t_{\frac{N}{2}}\right)+\mathrm{K}\right]-6 u\left(t_{\frac{N}{2}}\right)+ \\ & 4\left[u\left(t_{\frac{N}{2}}\right)-h_{\frac{N}{2}} u^{\prime}\left(t_{\frac{N}{2}}\right)+\frac{h_N^2}{2} u^{\prime \prime}\left(t_{\frac{N}{2}}\right)+\mathrm{K}\right] \\ & \left.-\left[\left(u\left(t_{\frac{N}{2}}\right)-h_{\frac{N}{2}} u^{\prime}\left(t_{\frac{N}{2}}\right)+\mathrm{K}\right)-h_{\frac{N}{2}-1}\left(u^{\prime}\left(t_{\frac{N}{2}}\right)-h_{\frac{N}{2}} u^{\prime \prime}\left(t_{\frac{N}{2}}\right)+\mathrm{K}\right)\right] \right\rvert\, \\ & \leq\left|2 h_1\right|\left|u^{\prime}\left(t_{\frac{N}{2}}\right)\right| \leq 2 h_1 C_1 \varepsilon^{-1 / 2} e_1(t)\end{aligned}$
Bounding exponential functions by constant, we have $\left|T_{i, u}\right| \leq C N^{-1} \beta^{-1 / 2} \ln N$, for $i=\frac{N}{2}$.
Hence proved.
This section validates, the effectiveness of the proposed method and it is implemented on two test problems. Maximum absolute error (MAE) is evaluated by using double mesh principle which is defined as:
$\begin{aligned} & E_{\mathcal{E}}^N=\max _{t_i \in \Psi_{\mathcal{E}}^N}\left|u_{\varepsilon}^N\left(t_i\right)-u_{\varepsilon}^{2 N}\left(t_i\right)\right| \text { and } E^N=\max E_{\varepsilon}^N .\end{aligned}$
where, $u_{\varepsilon}^N\left(t_i\right)$ is numerical solution at point $t_i$ for number of mesh points N and perturbation parameter $\varepsilon$ and $u_{\varepsilon}^{2 N}\left(t_i\right)$ is numerical solution at point $t_i$ for number of mesh points 2N.
The rate of convergence (ROC) is calculated through the given below formula:
$p^N=\log _2\left(\frac{E^N}{E^{2 N}}\right)$
Both the test problems have an interior and boundary layer. All the mathematical computations are performed in MATLAB software. The problems are solved for maximum absolute errors for various values of mesh points and perturbation parameter. Tables 1 and 2 show the maximum absolute error for example 1 and 2 respectively which is compared with example 1 for existing method [6] in Table 3 and it has been shown that the error obtained using current methodology is lesser than the method adopted in study [6]. Tables 4 and 5 represent the order of convergence for example 1 and 2 respectively. Figures 1 and 2 depict the of numerical solution of example 1 and 2 for mesh points N=64 and perturbation parameter $\varepsilon=2^{-6}$. These numerical results are in agreement with the theoretical results presented in this paper.
Example 1: Consider SPBVP (1-2) having:
$b(t)=1, g(t)=\left\{\begin{array}{rr}0.7, & \mathrm{t} \leq 0.5, \\ -0.6, & \mathrm{t}>0.5,\end{array} u(0)=1, u(1)=0\right.$
Figure 1. Graphical representation of numerical solution of example 1 for N=64 and $\varepsilon=2^{-6}$
Figure 2. Graphical representation of numerical solution of example 2 for N=64 and $\varepsilon=2^{-6}$
Table 1. MAE of example 1 for different values of N and $\varepsilon$
|
$\varepsilon / N$ |
64 |
128 |
256 |
512 |
1024 |
2048 |
|
$2^{-1}$ |
2.1500 E-03 |
1.0062 E-03 |
4.8534 E-04 |
2.3815 E-04 |
1.1794 E-04 |
5.8682 E-05 |
|
$2^{-2}$ |
3.7700 E-03 |
1.7569 E-03 |
8.4512 E-04 |
4.1406 E-04 |
2.0518 E-04 |
1.0220 E-04 |
|
$2^{-3}$ |
6.1148 E-03 |
2.8273 E-03 |
1.3689 E-03 |
6.7358 E-04 |
3.3409 E-04 |
1.6637 E-04 |
|
$2^{-4}$ |
9.0974 E-03 |
4.2614 E-03 |
2.0600 E-03 |
1.0125 E-03 |
5.0187 E-04 |
2.4985 E-04 |
|
$2^{-5}$ |
1.2967 E-02 |
6.0477 E-03 |
2.9167 E-03 |
1.4318 E-03 |
7.0929 E-04 |
3.5300 E-04 |
|
$2^{-6}$ |
1.7870 E-02 |
8.3162 E-03 |
4.0071 E-03 |
1.9662 E-03 |
9.7386 E-04 |
4.8462 E-04 |
|
$2^{-7}$ |
2.4221 E-02 |
1.1310 E-02 |
5.4624 E-03 |
2.6841 E-03 |
1.3304 E-03 |
6.6227 E-04 |
|
$2^{-8}$ |
3.2494 E-02 |
1.5345 E-02 |
7.4608 E-03 |
3.6790 E-03 |
1.8269 E-03 |
9.1030 E-04 |
|
$E^N$ |
3.2494 E-02 |
1.5345 E-02 |
7.4608 E-03 |
3.6790 E-03 |
1.8269 E-03 |
9.1030 E-04 |
Table 2. MAE of example 2 for different values of N and $\varepsilon$
|
$\varepsilon / N$ |
64 |
128 |
256 |
512 |
1024 |
2048 |
|
$2^0$ |
1.3462 E-03 |
6.6183 E-04 |
3.2959 E-04 |
1.6446 E-04 |
8.2142 E-05 |
4.1049 E-05 |
|
$2^{-1}$ |
2.3615 E-03 |
1.1489 E-03 |
5.7276 E-04 |
2.8593 E-04 |
1.4285 E-04 |
7.1395 E-05 |
|
$2^{-2}$ |
3.8059 E-03 |
1.8411 E-03 |
9.0801 E-04 |
4.5359 E-04 |
2.2669 E-04 |
1.1331 E-04 |
|
$2^{-3}$ |
5.5340 E-03 |
2.6563 E-03 |
1.3042 E-03 |
6.4622 E-04 |
3.2165 E-04 |
1.6060 E-04 |
|
$2^{-4}$ |
7.7134 E-03 |
3.7310 E-03 |
1.8345 E-03 |
9.0954 E-04 |
4.5286 E-04 |
2.2595 E-04 |
|
$2^{-5}$ |
1.0231 E-02 |
4.9529 E-03 |
2.4372 E-03 |
1.2090 E-03 |
6.0213 E-04 |
3.0047 E-04 |
|
$2^{-6}$ |
1.3170 E-02 |
6.3826 E-03 |
3.1441 E-03 |
1.5607 E-03 |
7.7754 E-04 |
3.8808 E-04 |
|
$2^{-7}$ |
1.6834 E-02 |
8.1928 E-03 |
4.0463 E-03 |
2.0113 E-03 |
1.0028 E-03 |
5.0069 E-04 |
|
$2^{-8}$ |
2.1672 E-02 |
1.0662 E-02 |
5.2942 E-03 |
2.6384 E-03 |
1.3171 E-03 |
6.5802 E-04 |
|
$E^N$ |
2.1672 E-02 |
1.0662 E-02 |
5.2942 E-03 |
2.6384 E-03 |
1.3171 E-03 |
6.5802 E-04 |
Table 3. MAE of example 1 by Farrell et al. [6] for different values of N and $\varepsilon$
|
$\varepsilon / N$ |
64 |
128 |
256 |
512 |
1024 |
2048 |
|
$2^{-1}$ |
4.3680 E-03 |
2.1820 E-03 |
1.0890 E-03 |
5.4200 E-04 |
2.6900 E-04 |
1.3200 E-04 |
|
$2^{-2}$ |
7.7270 E-03 |
3.8600 E-03 |
1.9260 E-03 |
9.5900 E-04 |
4.7600 E-04 |
2.3400 E-04 |
|
$2^{-3}$ |
1.2746 E-02 |
6.3680 E-03 |
3.1780 E-03 |
1.5830 E-03 |
7.8500 E-04 |
3.8600 E-04 |
|
$2^{-4}$ |
1.9557 E-02 |
9.7720 E-03 |
4.8760 E-03 |
2.4290 E-03 |
1.2050 E-03 |
5.9300 E-04 |
|
$2^{-5}$ |
2.8476 E-02 |
1.4233 E-02 |
7.1040 E-03 |
3.5380 E-03 |
1.7550 E-03 |
8.6400 E-04 |
|
$2^{-6}$ |
4.0484 E-02 |
2.0250 E-02 |
1.0109 E-02 |
5.0350 E-03 |
2.4980 E-03 |
1.2290 E-03 |
|
$2^{-7}$ |
5.7174 E-02 |
2.8642 E-02 |
1.4303 E-02 |
7.1250 E-03 |
3.5350 E-03 |
1.7390 E-03 |
|
$2^{-8}$ |
8.0544 E-02 |
4.0467 E-02 |
2.0223 E-02 |
1.0076 E-02 |
4.9990 E-03 |
2.4600 E-03 |
|
$E^N$ |
8.0544 E-02 |
4.0467 E-02 |
2.0223 E-02 |
1.0076 E-02 |
4.9990 E-03 |
2.4600 E-03 |
Table 4. ROC of example 1 for different values of N and $\varepsilon$
| $\varepsilon / N$ |
64 |
128 |
256 |
512 |
1024 |
|
$2^{-1}$ |
1.095 |
1.052 |
1.027 |
1.014 |
1.007 |
|
$2^{-2}$ |
1.102 |
1.056 |
1.029 |
1.013 |
1.006 |
|
$2^{-3}$ |
1.113 |
1.046 |
1.023 |
1.012 |
1.006 |
|
$2^{-4}$ |
1.094 |
1.049 |
1.025 |
1.012 |
1.006 |
|
$2^{-5}$ |
1.100 |
1.052 |
1.027 |
1.013 |
1.007 |
|
$2^{-6}$ |
1.104 |
1.053 |
1.027 |
1.014 |
1.007 |
|
$2^{-7}$ |
1.099 |
1.050 |
1.025 |
1.013 |
1.006 |
|
$2^{-8}$ |
1.082 |
1.040 |
1.020 |
1.010 |
1.005 |
Example 2: Consider SPBVP (1-2) having:
$\begin{aligned} & b(t)=\left\{\begin{array}{l}2 t+1, t \leq 0.5, \\ 3-2 t, t>0.5,\end{array}\right. \\ & g(t)=\left\{\begin{array}{l}-0.5, t \leq 0.5, \\ 0.5, t>0.5,\end{array}\right. \\ & u(0)=1, u(1)=g(0) .\end{aligned}$
Table 5. ROC of example 2 for different values of N and $\varepsilon$
|
$\varepsilon / N$ |
64 |
128 |
256 |
512 |
1024 |
|
$2^{-1}$ |
1.039 |
1.004 |
1.002 |
1.001 |
1.001 |
|
$2^{-2}$ |
1.048 |
1.020 |
1.001 |
1.001 |
1.000 |
|
$2^{-3}$ |
1.059 |
1.026 |
1.013 |
1.007 |
1.002 |
|
$2^{-4}$ |
1.048 |
1.024 |
1.012 |
1.006 |
1.003 |
|
$2^{-5}$ |
1.047 |
1.023 |
1.011 |
1.006 |
1.003 |
|
$2^{-6}$ |
1.045 |
1.022 |
1.010 |
1.005 |
1.003 |
|
$2^{-7}$ |
1.039 |
1.018 |
1.008 |
1.004 |
1.002 |
|
$2^{-8}$ |
1.023 |
1.010 |
1.005 |
1.002 |
1.001 |
An effective hybrid numerical method is proposed to solve SPBVP of reaction-diffusion type with DST. The CSM is applied on fine regions and a classical central FDM is applied on coarse regions based on Shishkin mesh. The present study improves the result by reducing the MAE and increasing ROC. Proposed method is applied on two test problems and numerical results are compared with the exiting method [6]. It gives first order convergence result which satisfy the theoretical results. In future, researchers can use some other numerical methods like higher order spline method to reduce the MAE and increase the ROC.
|
$\mathcal{E}$ |
Perturbation parameter |
|
$w$ |
Point of discontinuity |
|
$C^1$ |
Continuously first differentiable function |
|
$C^2$ |
Twice Continuously differentiable function |
|
$L_{\varepsilon}$ |
Differential operator |
|
$C$ |
Constant |
|
$\mu_1$ and $\mu_2$ |
Transition parameter |
|
$N$ |
Number of mesh points |
|
$p^N$ |
Rate of convergence |
[1] Miller, J.J.H., O'Riordan, E., Shishkin, G.I., Wang, S. (2000). A parameter-uniform schwarz method for a singularly perturbed reaction-diffusion problem with an interior layer. Applied Numerical Mathematics, 35(4): 323-337. https://doi.org/10.1016/S0168-9274(99)00140-3
[2] Farrell, P., Hegarty, A., Miller, J.M., O'Riordan, E., Shishkin, G.I. (2000). Robust Computational Techniques for Boundary Layers. CRC Press.
[3] Miller, J.J.H., O'Riordan, E., Shishkin, G.I. (2012). Fitted Numerical Methods for Singular Perturbation Problems. World-Scientific, Singapore.
[4] Roos, H.G., Stynes, M., Tobiska, L. (1996). Numerical methods for singularly perturbed differential equations. Springer Series in Computational Mathematics (SSCM). https://doi.org/10.1007/978-3-662-03206-0
[5] O'Malley, R.E. (1991). Singular perturbation methods for ordinary differential equations. In Applied Mathematical Sciences (AMS), New York: Springer, p. 225. https://doi.org/10.1007/978-1-4612-0977-5
[6] Farrell, P.A., Miller, J.J.H., O'Riordan, E., Shishkin, G.I. (2000). Singularly perturbed differential equations with discontinuous source terms. Analytical and Numerical Methods for Convection-dominated and Singularly Perturbed Problems, 23-32.
[7] Roos, H.G., Zarin, H. (2002). A second-order scheme for singularly perturbed differential equations with discontinuous source term. Journal of Numerical Mathematics, 10(4): 275-289. https://doi.org/10.1515/JNMA.2002.275
[8] Shanthi, V., Ramanujam, N., Natesan, S. (2006). Fitted mesh method for singularly perturbed reaction-convection-diffusion problems with boundary and interior layers. Journal of Applied Mathematics and Computing, 22: 49-65. https://doi.org/10.1007/BF02896460
[9] Chandru, M., Shanthi, V. (2014). A boundary value technique for singularly perturbed boundary value problem of reaction-diffusion with non-smooth data. In Journal of Engineering Science and Technology, 32-45.
[10] Chandru, M., Prabha, T., Shanthi, V. (2015). A hybrid difference scheme for a second-order singularly perturbed reaction-diffusion problem with non-smooth data. International Journal of Applied and Computational Mathematics, 1: 87-100. https://doi.org/10.1007/s40819-014-0004-8
[11] Chandru, M., Shanthi, V. (2015). Fitted mesh method for singularly perturbed robin type boundary value problem with discontinuous source term. International Journal of Applied and Computational Mathematics, 1: 491-501. https://doi.org/10.1007/s40819-015-0030-1
[12] De Falco, C., O'Riordan, E. (2010). Interior layers in a reaction-diffusion equation with a discontinuous diffusion coefficient. International Journal of Numerical Analysis and Modeling, 7(3): 444-461.
[13] Cen, Z.D., Le, A.B., Xu, A.M. (2017). A high-order finite difference scheme for a singularly perturbed reaction-diffusion problem with an interior layer. Advances in Difference Equations, 2017(1): 1-14. https://doi.org/10.1186/s13662-017-1268-1
[14] Natesan, S., Bawa, R.K. (2007). Second-order numerical scheme for singularly perturbed reaction-diffusion robin problems. Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM), 2(3-4): 177-192.
[15] Natesan, S., Deb, B.S. (2007). A robust computational method for singularly perturbed coupled system of reaction-diffusion boundary-value problems. Applied Mathematics and Computation, 188(1): 353-364. https://doi.org/10.1016/j.amc.2006.09.120
[16] Bawa, R.K., Natesan, S. (2005). A computational method for self-adjoint singular perturbation problems using quintic spline. Computers & Mathematics with Applications, 50(8-9): 1371-1382. https://doi.org/10.1016/j.camwa.2005.04.017
[17] Kadalbajoo, M.K., Bawa, R.K. (1996). Variable-mesh difference scheme for singularly-perturbed boundary-value problems using splines. Journal of Optimization Theory and Applications, 90(2): 405-416. https://doi.org/10.1007/BF02190005
[18] Kadalbajoo, M.K., Bawa, R.K. (1993). Third-order variable-mesh cubic spline methods for singularly-perturbed boundary-value problems. Applied Mathematics and Computation, 59(2-3): 117-129. https://doi.org/10.1016/0096-3003(93)90084-R
[19] Natesan, S., Vigo-Aguiar, J.A., Ramanujam, N. (2003). A numerical algorithm for singular perturbation problems exhibiting weak boundary layers. Computers & Mathematics with Applications, 45(1-3): 469-479. https://doi.org/10.1016/S0898-1221(03)80031-7
[20] Kadalbajoo, M., Patidar, K. (2003). Variable mesh spline in compression for the numerical solution of singular perturbation problems. International Journal of Computer Mathematics, 80(1): 83-93. https://doi.org/10.1080/00207160304656
[21] Chawla, M.M., Subramanian, R. (1987). A new fourth order cubic spline method for non-linear two-point boundary value problems. International Journal of Computer Mathematics, 22(3-4): 321-341. https://doi.org/10.1080/00207168708803601
[22] Kadalbajoo, M.K., Aggarwal, V.K. (2005). Fitted mesh B-spline collocation method for solving self-adjoint singularly perturbed boundary value problems. Applied Mathematics and Computation, 161(3): 973-987. https://doi.org/10.1016/j.amc.2003.12.078
[23] Alberg, J.H., Nilson, E.N., Walsh, J.L. (1967). The theory of splines and their applications. New York: Academic.
[24] Prenter, P.M. (2008). Splines and variational methods. Courier Corporation.
[25] Micula, G., Micula, S. (1999). Handbook of Splines. Springer Science & Business Media.
[26] Tamilselvan, A., Ramanujam, N., Shanthi, V. (2007). A numerical method for singularly perturbed weakly coupled system of two second order ordinary differential equations with discontinuous source term. Journal of Computational and Applied Mathematics, 202(2): 203-216. https://doi.org/10.1016/j.cam.2006.02.025
