In this paper, we consider the following singularly perturbed problem  ;
subject to the following boundary conditions,
where and are two small perturbation parameters; and are sufficiently smooth functions for ; a, b, α, and β are real constants. In general, the solution may exhibit two boundary layers of exponential type at both end points .
Different applications in science and engineering consider these kinds of problems that describe complicated physical and chemical models such as heat transfer problems, Navier-Stokes flow with large Reynolds numbers, chemical reactor theory, convection-diffusion processes, geophysics, aerodynamics, reaction-diffusion processes, quantum mechanics and optimal control, etc. Its solution exhibits two layers at the two endpoints of the domain. The nature of the two-parameter problem was asymptotically examined by  . It was found that layer-adapted meshes have been required to obtain a uniformly convergent method no matter how small the perturbation parameter see  for more details.
Many numerical methods have been developed for the solution of two layer boundary value problems, such as described in  ,  ,  and  for one parameter singularly perturbed boundary value problems and with two small parameters are considered in  ,  ,  and  , but on a Shishkin-type mesh. Vulanovic  considered Shishkin and Bakhvalov meshes but assumed
with p > 0. Dag and Sahin presented a numerical solution of singularly
perturbed boundary value problems, using finite element method  . Their collocation method was applied with quadratic and cubic B-spline base functions over the geometrically graded mesh of the solution domain. In 2010, Rashidina and Mohammadi  considered the self-adjoint singularly perturbed two-point boundary value problems. Ramadan et al. (2007)  developed quintic nonpolynomial spline methods for the numerical solution of fourth order two-point boundary value problems. A second order monotone numerical method was constructed by Gracia et al. in 2006  for a singularly perturbed ordinary differential equation with two small parameters affecting the convection and diffusion terms. The monotone operator was combined with a piecewise uniform Shishkin mesh. Kadalbajoo and Yadaw  presented a B-spline collocation method for solving a class of two-parameter singularly perturbed boundary value problems. They used B-spline collocation method on piecewise-uniform Shishkin mesh, which leads to a tridiagonal linear system. Their method was shown to have a uniform convergence of second order.
For a more reasonable analysis of the system response, uncertainty should be involved. When the exact value of a quantity is unknown, its approximation and corresponding degree of uncertainty can be conveyed via an interval, which estimates a range of possible values expected to include, firstly suggested by several mathematicians for bounding round-off errors, the interval analysis is fully developed by Moore  . We are going to use MCS as a validating tool for the proposed method. MCS is the simplest method for treating any randomness in a system  . The method basically depends on generating a set of realization of the random parameter, then a unique solution is defined by carrying out the deterministic solver for each of these realizations. In 1967, Stein  generalized his model which incorporates stochastic effects due to neuronal excitations to handle a distribution of post-synaptic potential amplitudes and used the Monte-Carlo technique for approximating the solution. David Edwards  developed a multi-region FDM technique for a particular singularly perturbed boundary value problem and this method was based on Monte Carlo techniques.
The main contribution of this paper is to develop a new spline method based on a Shishkin mesh discretization for obtaining an approximation for the solution of two-layer boundary value problems. As the perturbation parameter is not deterministic, therefore the interval analysis is considered to estimate the solution range. The validation of the developed solver will be done by comparing the exact and the approximate solutions of the proposed method and the MCS results. The convergence analysis also is presented numerically and shows that the presented method is almost second-order.
The paper is organized as follows: In Section 2, we derive our spline scheme. Mesh strategy based on a Shishkin mesh is presented in Section 3. In Section 4, we present interval analysis and sensitivity measures. Numerical results are discussed in Section 5. Section 6 is devoted to the final conclusions, while future work is provided in section 7. Finally, Section 8 is dedicated to references.
2. Derivation of Exponential Spline
We discretize the solution region such that . Where N is the number of mesh points. Let be the mesh size and the mesh ratio
. When the mesh reduces to a uniform mesh, . The interpolating exponential spline approximation function can be defined as  :
where , , are constants and τ is a free parameter such that the non-polynomial spline (3) reduces to usual cubic spline when τ approaches to zero  , which satisfies the following conditions:
2) . (4)
The algebraic manipulations of Equations (3) and (4) yield the following expressions:
4) . (5)
From the aspect of the first derivative continuity at the mesh points yields the expression for the determination of where . We can get the following exponential spline identity relation:
Note that, the exponential spline relation (6) is consistent with the standard variable-mesh cubic spline if , hence  . If we choose , this method is known as the Numerov method  .
3. Mesh Selection Strategy
We consider the simplest possible non-uniform mesh, namely a piecewise-uniform mesh proposed by Shishkin  . The domain is into three sub-domains
where the transition parameters are given by:
where and are two solutions of the characteristic equation:
The quantity describes the boundary layer at , while characterizes the layer at , and
We take N/4, N/2 and N/4 mesh points, respectively in and . Denote the step sizes in each subinterval by , and , respectively. Accordingly, the resulting piecewise-uniform Shishkin mesh is represented by:
4. Interval and Sensitivity Analysis
Since uncertainty will be considered, the interval analysis can be used as descriptive measures of uncertainty in quantitative values. Hence, the perturbation parameter is not deterministic, the solution has to be defined as a range based on the interval of the parameter. Therefore, the upper and lower bounds of the perturbation parameter can be written as:
where is the upper value, is the lower value and is the central value. Then the fluctuation range of solution could be estimated.
Sensitivity measures can be conducted using different techniques for example One-at-a-Time Sensitivity Measures (±SD), the Sensitivity Index (SI), the Importance Index (II), Differential Sensitivity Analysis (PD), etc. We estimated the sensitivity measures using the following methods One-at-a-Time Sensitivity Measures (±SD), the Sensitivity Index, and the Differential Sensitivity Analysis. One-at-a-Time Sensitivity Measures (±SD) is considered the simplest―but powerful―method for conducting sensitivity analysis, which estimates the variation of the solution as the perturbation parameter is increased by a factor of its standard deviation or in other words a percentage of its mean value  .
The Sensitivity Index (SI), another simple method of estimating the sensitivity measure is to calculate the relative solution difference when varying one input parameter from its minimum value to its maximum value, which provides a good indication of parameter and model variability. The SI is calculated using, , where Ymin and Ymax represent the minimum and maximum solution values, respectively, resulting from varying the perturbation parameter over its entire range  .
Differential Sensitivity Analysis (PD) method considers all random parameters equal to their mean values and partial differentiation of the system with respect to the random parameters should be done. The sensitivity coefficient for a specific parameter ?perturbation parameter can be measured from the partial derivative relation . The results are normalized by multiplying the derivatives by the ratio of the parameter value to the solution for the mean value, sensitivity coefficient  .
5. Numerical Example
We consider the following reaction-diffusion problem; see   and  :
whose exact solution is given by:
The estimated maximum error and the rate of convergence are computed by the formulas:
Table 1 shows the maximum absolute error and the order of convergence for various values of the perturbation parameter ℇ. The results obtained using the current method are very accurate compared with the analytical solution and gives the order of convergence 2 even for small values of ℇ, and it is shown in Figure 1 that the exact and the approximate solutions are very close. Furthermore, since the problem is singularly perturbed its solution possesses layers along the boundary of the domain which is occurred in the form of sharp boundary-layers in Figure 1 at .
Figures 2-4 show the comparison of the solution between the proposed method and the MCS (50000 samples), as the perturbation parameter changes by 20%.
Table 1. Maximum absolute errors and the order of convergence for .
Figure 1. Deterministic Case: exact and approximate solutions for , .
Figure 2. Random Case: Mean of the solution of the proposed method and the MCS results, , .
Figure 3. Random Case: Upper limit of the solution range of the proposed method and the MCS results, , .
Figure 4. Random Case: limit of the solution range of the proposed method and the MCS results, , .
The upper, centre and lower solutions of the two methods are very close which shows the accuracy and validation of the proposed method. Further sensitivity measures have been conducted such as the SI method gives an index of value 0.0039%, which is very small. The differential method indicates that the sensitivity coefficient is also very small by value 1.34E−15 at the mid-point of the scheme. Therefore, the solution is not sensitive to the changes in the perturbation parameter.
A numerical method based on exponential spline with Shishkin mesh discretization is combined with interval analysis perspective to evaluate the range of the solution for the singularly perturbed two-point boundary value problems with uncertain parameter. The numerical results show that the present method approximates the solution very well compared with the exact solution. Therefore, the proposed method is almost second-order uniformly convergent with respect to the perturbation parameter. MCS are used to prove the validation and the accuracy of the proposed method. Sensitivity analysis has been conducted using different methods and it is found that the solution is not sensitive to the perturbation parameter.
7. Future Work
This paper proposed a new numerical spline method combined with interval analysis to solve singularly perturbed boundary value problems. We used interval analysis to estimate the solution range as the perturbation parameter is not deterministic and compared the results with the MCS method to validate the proposed method. In the future, we shall use another stochastic method to deal with the uncertainty, which appears in the perturbation parameter for example polynomial chaos expansion.
 Zahra, W.K. and Van Daele, M. (2015) Uniformly Convergent Discrete Spline Scheme on a Shishkin Mesh for the Singular Perturbation Boundary Value Problem. Proceedings of the 15th International Conference on Mathematical Methods in Science and Engineering, CMMSE 2015, 6-10 July 2015, Cádiz, 1261-1268.
 O’Malley, R.E. (1967) Singular Perturbations of Boundary Value Problems for Linear Ordinary Differential Equations Involving Two Parameters. Journal of Mathematical Analysis and Applications, 19, 291-308.
 Zahra, W.K. and Van Daele, M. (2018) Discrete Spline Solution of Singularly Perturbed Problem with Two Small Parameters on a Shishkin-Type Mesh. Computational Mathematics and Modeling Journal. (To be published)
 Kumar, D., Yadaw, A.S. and Kadalbajoo, M.K. (2013) A Parameter-Uniform Method for Two Parameters Singularly Perturbed Boundary Value Problems via Asymptotic Expansion. Applied Mathematics & Information Sciences, 7, 1525-1532.
 Duvnjakovic, E., Karasuljic, S., Pasic, V. and Zarin, H. (2014) A Uniformly Convergent Difference Scheme on a Modified Shishkin Mesh for the Singular Perturbation Boundary Value Problem. arXiv Prepr. arXiv1411.4323
 Ramadan, M.A., Lashien, I.F. and Zahra, W.K. (2007) The Numerical Solution of Singularly Perturbed Boundary Value Problems Using Nonpolynomial Spline. International Journal of Pure and Applied Mathematics, 41, 883-896.
 Zahra, W.K., El-Azab, M.S. and El Mhlawy, A.M. (2014) Spline Difference Scheme for Two-Parameter Singularly Perturbed Partial Differential Equations. Journal of Applied Mathematics & Informatics, 32, 185-201.
 Natesan, S., Gracia, J.L. and Clavero, C. (2004) Singularly Perturbed Boundary-Value Problems with Two Small Parameters: A Defect Correction Approach. Proceedings of the International Conference on Boundary and Interior Layers, Computational and Asymptotic Methods, BAIL, Bail, 1-6.
 Kumar, M. and Rao, S.C.S. (2010) High Order Parameter-Robust Numerical Method for Singularly Perturbed Reaction-Diffusion Problems. Applied Mathematics and Computation, 216, 1036-1046.
 Rashidinia, J. and Mohammadi, R. (2010) Non-Polynomial Spline Approximations for the Solution of Singularlyperturbed Boundary Value Problems. TWMS Journal of Pure and Applied Mathematics, 1, 236-251.
 Ramadan, M.A., Lashien, I.F. and Zahra, W.K. (2009) Quintic Nonpolynomial Spline Solutions for Fourth Order Two-Point Boundary Value Problem. Communications in Nonlinear Science and Numerical Simulation, 14, 1105-1114.
 Gracia, J.L., O’Riordan, E. and Pickett, M.L. (2006) A Parameter Robust Second Order Numerical Method for a Singularly Perturbed Two-Parameter Problem. Applied Numerical Mathematics, 56, 962-980.
 Kadalbajoo, M.K. and Yadaw, A.S. (2008) B-Spline Collocation Method for a Two-Parameter Singularly Perturbed Convection-Diffusion Boundary Value Problems. Applied Mathematics and Computation, 201, 504-513.
 Edwards, D. (2009) High Precision Calculations of One Dimension Singularly Perturbed Boundary Value Problems Using Multi Region FDM. Proceedings of the International MultiConference of Engineers and Computer Scientists, Vol. II, Hong Kong, 18-20 March 2009.
 Zahra, W.K. (2011) Finite-Difference Technique Based on Exponential Splines for the Solution of Obstacle Problems. International Journal of Computer Mathematics, 88, 3046-3060.
 Tirmizi, I.A., Fazal-i-Haq and Siraj-ul-Islam (2008) Non-Polynomial Spline Solution of Singularly Perturbed Boundary-Value Problems. Applied Mathematics and Computation, 196, 6-16.
 Mohanty, R.K., Nayak, S. and Khan, A. (2017) Non-Polynomial Cubic Spline Discretization for System of Non-Linear Singular Boundary Value Problems Using Variable Mesh. Advances in Difference Equations, 2017, 327.
 Henrici, P. (1962) Discrete Variable Methods in Ordinary Differential Equations. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschriftfür Angewandte Mathematik und Mechanik.
 Shishkin, G.I. (1988) Grid Approximation of Singularly Perturbed Parabolic Equations with Internal Layers. Russian Journal of Numerical Analysis and Mathematical Modelling, 3, 393-408.
 Hoffman, F.O. and Miller, C.W. (1983) Uncertainties in Environmental Radiological Assessment Models and Their Implications. Annual Meeting of the National Council of Radiation Protection and Measurements, Washington DC, 6 April 1983.
 Khan, A. and Khandelwal, P. (2014) Non-Polynomial Sextic Spline Solution of Singularly Perturbed Boundary-Value Problems. International Journal of Computer Mathematics, 91, 1122-1135.
 Zahra, W.K. and El Mhlawy, A.M. (2013) Numerical Solution of Two-Parameter Singularly Perturbed Boundary Value Problems via Exponential Spline. Journal of King Saud University-Science, 25, 201-208.