Fourth-order fitted mesh scheme for semilinear singularly perturbed reaction–diffusion problems

Abstract

Objective

The main purpose of this work is to present a fourth-order fitted mesh scheme for solving the semilinear singularly perturbed reaction–diffusion problem to produce more accurate solutions.

Results

Quasilinearization technique is used to linearize the semilinear term. The scheme is formulated with discretizing the solution domain piecewise uniformly and then replacing the differential equation by finite difference approximations. This gives the system of difference algebraic equations and is solved by the Thomas algorithm. Convergence analysis are investigated using solution bound and the truncation error bound. Numerical illustrations are investigated to support the theoretical results and the method’s applicability. The method produces a more accurate solution than some existing methods in the literature.

Introduction

Any differential equation in which the coefficient of highest order derivative is a small positive parameter together with the prescribed conditions is called singularly perturbed problem, [1]. This problem arise very frequently in diversified fields of applied mathematics and engineering; for instance fluid mechanics, elasticity, hydrodynamics, quantum mechanics, chemical-reactor theory, aerodynamics, plasma dynamics, modeling of semiconductor devices, diffraction theory and reaction–diffusion processes and many other allied areas [2,3,4,5,6,7,8]. Hence, due to the importance of these problems in real life situations, the need to develop numerical methods for approximation its solution is advantageous.

Singularly perturbed problems broadly categorized into reaction–diffusion and convection–diffusion types, [9,10,11,12,13,14,15, 22]. Thes can be further classified according to the type of layers (boundary and/or interior layers), location (left or right or twin), type of defined conditions like initial or boundary (Dirichlet, robin or mixed conditions). More particularly, the singularly perturbed reaction–diffusion boundary valued problems are categorized into linear and nonlinear problems exhibiting twin (both left and right) boundary layers. These types of problems occur frequently in fluid mechanics and other several fields of applied mathematics, [16,17,18,19,20,21,22].

As both books in [10, 15], explain, the region at which the solution of singularly perturbed problems change very quickly in certain small regions called layer region. It is well known that when the perturbation parameter is small enough, classical numerical methods fails to produce the required accurate solution for these problems. Thus, it should be important to develop appropriate numerical methods for such problems. There are several numerical methods suggested by various authors for solving the families of the linear singularly perturbation reaction–diffusion problems. Some and most of them are higher order (fourth, sixth, eighth, and tenth order) compact finite difference methods to solve different types singularly perturbed reaction–diffusion problems [2,3,4,5]. Most authors in these literatures developed the methods by restricting the criteria for the constant coefficients of the reaction term within the problem. Further, scholars in [2, 13], have presented fourth and sixth-order stable central difference method, respectively, for solving singularly perturbed two-point boundary value problem. This type of problem shares the basic behavior of singularly perturbed reaction–diffusion boundary value problem. Authors in [23] presented a numerical approach to solving singularly perturbed semilinear convection–diffusion problems. The nonlinear part of the problem is linearized via the quasilinearization technique.

Recently, a numerical scheme for solving the semilinear singularly perturbed reaction–diffusion problems and a numerical technique for solving a class of nonlinear singularly perturbed boundary value problems have been presented respectively, in [2, 14]. They have introduced a basic and computational approach scheme based on Numerov’s type on uniform mesh. They indicated that the method is uniformly convergence, according to the discrete maximum norm, independently of the perturbation parameter. Also, authors in [3], deliberate a numerical scheme based on Sinc collocation method to solve a class of nonlinear singularly perturbed boundary value problems which exhibit twin boundary layers. The Sinc method can control the oscillations in computed solutions at boundary layer regions naturally because the distribution of Sinc points is denser at near the boundaries. The numerical results support the theoretical results and illustrate the validity and accuracy of the method compared with the results in the existing methods. So far different finite difference methods have been adopted for solving the semilinear singularly perturbed reaction–diffusion problems. However, the obtained accurate solution and the existing rate of convergence are attracts remarkable attention to be improve. Thus, the main drawbacks to most of the presented methods are using uniform mesh of the solution domain, restricting the coefficient reaction term to constant function only, the method gives accurate solution when the mesh size of the solution domain and the values of the perturbation parameter are comparable. Moreover, most presented numerical schemes for solving the considered problem are limited to the second order convergent.

Therefore, in this paper, the main objective is to present a fourth-order fitted mesh scheme that works for variable coefficient of the reaction term for solving the semilinear singularly perturbed boundary value problems. Accordingly, in this effort, we have been formulated, analyze and implement the fourth-order fitted mesh scheme that produces a more accurate solution for solving semilinear singularly perturbed reaction–diffusion problems.

Main text

Description of the method

This paper deals with singularly perturbed semilinear reaction–diffusion problem:

$$ - \varepsilon y^{\prime\prime}(x) + q(x)y(x) + g(x,y(x)) = 0,\quad \forall x \in (0,1), $$

(1)

subject to the conditions

$$ y(0) = A,\quad y(1) = B, $$

(2)

where $\varepsilon ,\,\,\,0 < \varepsilon < < 1$ is perturbation parameter, $A$ and $B$ are given constant numbers. Assume that the functions $q(x)$ and the nonlinear term $g(x,y(x))$ is given sufficiently smooth with $g(x,y(x)) \in C[0,1],$ [0,1], and

$$ q(x) + \frac{\partial g(x,y(x))}{{\partial y}} \ge \beta > 0,\quad \forall x \in [0,1], $$

(3)

for some constant $\beta > 0,$ to ensure the existence and unique solution with dual boundary layers near the two end points of the solution domain [14]. The reduced problem of Eq. (1) is obtained by setting $\varepsilon = 0,$ gives:

$$ q(x)y(x) + g(x,y(x)) = 0\,. $$

(4)

With these conditions Eqs. (1) and (2) and the reduced problem in Eq. (4) have a unique solution. The unique solution to Eqs. (1) and (2) exhibits two boundary layers at the end of the interval $\overline{\Omega } = [0,\,1]$, as the perturbation parameter $\varepsilon$ approaches to zero [2, 20, 22].

Quasilinearization technique used to transform the semilinear singularly perturbed reaction–diffusion problem into a sequence of linear equations, [14]. We choose a reliable initial approximation for the function $y^{(0)} (x)$ in $g(x,y(x))$ as:

$$ y^{(0)} (x) = mx + b,\, $$

(5)

where $m$ and $b$ are arbitrary constants determined using Eq. (2). Hence, Eq. (5) becomes:

$$ y^{(0)} (x) = (B - A)x + A. $$

(6)

By Taylor series expansion, we expand $g(x,y(x))$ around the chosen initial approximation:

$$ g(x,y^{(k + 1)} (x)) \simeq g(x,y^{(k)} (x)) + (y^{(k + 1)} - y^{k} )\frac{{\partial g^{(k)} }}{\partial y}\,|_{{(x,y^{(k)} (x))}} + \,\,...\,, $$

(7)

where $k = 0,1,2,...$ is the number of iterations index. Substituting Eq. (7) into Eq. (1) and (2), we have:

$$ \begin{aligned} - \varepsilon y^{\prime \prime (k + 1)} (x) & + q(x)y^{(k + 1)} (x) + \frac{{\partial g^{(k)} }}{\partial y}|_{{(x,y^{(k)} (x)}} y^{(k + 1)} (x) \hfill \\ & = - g(x,y^{(k)} (x)) + y^{(k)} (x)\frac{{\partial g^{(k)} }}{\partial y}|_{{_{{(x,y^{(k)} (x)}} }} + \cdots \,\,, \hfill \\ \end{aligned} $$

(8)

subject to the boundary conditions

$$ y^{(k + 1)} (0) = A,\,\,\,\,\,\,\,y^{(k + 1)} (1) = B\, $$

(9)

From now, Eq. (8) is linear in $y^{(k + 1)} (x).$ Thus, we solve the sequence of Eq. (8) in place of the semilinear problem in Eq. (1) by fourth-order fitted mesh scheme, which will be introduced in the next section. For the solution of the semilinear boundary value problem, we require that

$$ \mathop {\max \,}\limits_{k \to \infty } y^{(k)} (x) = y^{(*)} (x),\,\, $$

(10)

where $y^{(*)} (x)$ is the solution of the semilinear problem. Numerically, we require that:

$$ |y^{(k + 1)} (x) - y^{(k)} (x)|\,\, < \,\lambda ,\, $$

(11)

where $\lambda$ is a small tolerance chosen by us. Then $y^{(k + 1)} (x)$ is the approximate solution of the semilinear problem.

Properties of continuous problem

For the sake of simplicity, at first iteration or $(k = 0)$, let us denote Eqs. (8) and (9) in the form of:

$$ Lu(x) = - \varepsilon u^{\prime\prime}(x) + p(x)u(x) = f(x),\,\,\,x \in \Omega : = (0,1)\,, $$

(12)

where the coefficient of reaction term and the source terms are respectively:

$p(x) = q(x) + \frac{\partial g(x,y(x))}{{\partial y}}|_{{(x,y^{(0)} (x))}} \, \ge \,\,\beta > 0,$ and $f(x) = - g(x,y^{(0)} (x)) + y^{(0)} (x)\frac{\partial g(x,y(x))}{{\partial y}}|_{{(x,y^{(0)} (x))}} ,$and $y^{(k + 1)} (x) = u(x)$ with the operator $L$, such that:

$$ u(0) = A,\,\,\,\,\,\,u(1) = B\,. $$

(13)

We present some important properties for the solution of Eqs. (12) and (13) which will be useful in the subsequent section for the analysis of relevant numerical solutions.

Lemma 1

(Continuous Maximum Principle), [14], Assume that $v(x)$ is sufficiently smooth function which satisfies $v(0) \ge 0$ and $v(1) \ge 0$. Then $Lv(x) \ge 0,\,\,\,\,0 < x < 1,$ implies that $v(x) \ge 0$ $\forall x \in [0,1].$

Proof

Let $v$ be a value such that $v(x^{*} ) = \mathop {\min }\limits_{x \in \Omega } v(x),$ and assume that $v(x^{*} ) < 0.$

Clearly, $x^{*} \notin \{ 0,1\} ,$ and therefore, $v^{\prime}(x^{*} ) = 0$ and $v^{\prime\prime}(x) \ge 0.$ Moreover, there is

$$ Lv(x^{*} ) = - \varepsilon v^{\prime\prime}(x^{*} ) + p(x)v(x) < 0 $$

which is a contradiction. It follows that $v(x^{*} ) \ge 0$ and thus, $v(x) \ge 0,\,\,\,\forall x \in \Omega$.$\square$

Lemma 2

(Uniform Stability Estimate), [14], Let $u(x)$ be the solution of Eqs. (4.12) and (4.13). Then, we have:

$$ ||u(x)||\,\, \le \,\,\beta^{ - 1} ||f||\, + \max (|A|,|B|),\,\,\,\forall x \in [0,1]. $$

Proof

We construct two barrier functions $\psi^{ \pm }$ defined by.

$$ \psi^{ \pm } = \beta^{ - 1} ||f|| + \max (|A|,|B|) \pm u(x) $$

Then it can be said that

$$ \begin{gathered} \psi^{ \pm } (0) = \beta^{ - 1} ||f|| + \max (|A|,\,|B|) \pm u(0) = \beta^{ - 1} ||f|| + \max (|A|,\,|B|) \pm A \ge 0; \hfill \\ \psi^{ \pm } (1) = \beta^{ - 1} ||f|| + \max (|A|,\,|B|) \pm u(1)\, = \beta^{ - 1} ||f|| + \max (|A|,\,|B|) \pm B\, \ge 0; \hfill \\ \end{gathered} $$

It follows that

$$ \begin{aligned} L\psi^{ \pm } (x) & = - \varepsilon (\psi^{ \pm } (x))^{\prime\prime} + p(x)\psi^{ \pm } (x) \hfill \\ & = - \varepsilon [\beta^{ - 1} ||f|| + \max (|A|,|B|) \pm u(x)]^{\prime\prime} + p(x)[\beta^{ - 1} ||f|| + \max (|A|,|B|) \pm u(x)] \hfill \\ & = \pm ( - \varepsilon u^{\prime\prime}(x)) \pm p(x)u(x) + p(x)[\beta^{ - 1} ||f|| + \max (|A|,|B|)] \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{aligned} $$

Because, $[\beta^{ - 1} ||f|| + \max (|A|,|B|)]^{\prime\prime} = 0$

$$ \begin{aligned} L\psi^{ \pm } (x) & = \pm [ - \varepsilon u^{\prime\prime}(x)) + p(x)u(x)] + p(x)[\beta^{ - 1} ||f|| + \max (|A|,|B|)] \hfill \\ & = \pm Lu(x) + p(x)[\beta^{ - 1} ||f|| + \max (|A|,|B|)] \hfill \\ & = p(x)[\beta^{ - 1} ||f|| + \max (|A|,|B|)] \pm Lu(x) \hfill \\ \end{aligned} $$

Since, $p(x) \ge \beta > 0,\,\,\,\,\,||f||\, \ge \,f(x),\,\,\,\,{\text{and}}\,\,\,\,Lu(x) = f(x)$

$$ L\psi^{ \pm } (x) = p(x)[\beta^{ - 1} ||f|| + \max (|A|,|B|)] \pm f(x)\,\,\, \ge \,\,0. $$

Therefore, $L\psi^{ \pm } (x) \ge 0$. Thus, using Lemma 1, we get $\psi^{ \pm } (x) \ge 0,\,\,\,\,\forall x \in [0,1]$ this gives the required estimate. Further, we derive the bounds for the solution $u(x)$ and its derivative by the following estimate. $\square$

Lemma 3

(Boundedness), [14], Let $u(x)$ be the solution of Eqs. (4.12) and (4.13), with $p(x)$ and $f(x)$ are given sufficiently smooth functions. Then the following estimates hold.

$$ ||u(x)||_{\infty } \le C,\,\,\,\,\,0 \le x \le 1 $$

(14)

$$ |u^{\prime}(x)|\,\, \le C\left\{ {1 + \frac{1}{\sqrt \varepsilon }\left( {e^{{ - \sqrt {\frac{\alpha }{\varepsilon }} x}} + e^{{ - \sqrt {\frac{\alpha }{\varepsilon }} (1 - x)}} } \right)} \right\}\, $$

(15)

Proof

Applying Lemma 2, to Eqs. (12) and (13), we have Eq. (14).

$$ Lv(x) = \phi (x). $$

(16)

$$ v(0) = O\left( {\frac{1}{\sqrt \varepsilon }} \right) = v(1), $$

(17)

where $v(x) = u^{\prime}(x).$

$$ \phi (x) = f^{\prime}(x) - p^{\prime}(x)u(x). $$

(18)

The solution of the problem in Eq. (16) and (17) has the following forms:

$$ v(x) = v_{0} (x) + v_{1} (x) $$

(19)

where the functions $v_{0} (x)$ and $v_{1} (x)$ are the solutions of the following problems respectively

$$ \left\{ \begin{gathered} Lv_{0} (x) = \phi (x),\,\,\,\,\,\,0 < x < 1, \hfill \\ v_{0} (0) = v_{1} (1) = 0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered} \right. $$

(20)

and

$$ \left\{ \begin{gathered} Lv_{1} (x) = 0,\,\,\,0 < x < 1, \hfill \\ v_{0} (0) = v_{1} (1) = 0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \hfill \\ \end{gathered} \right. $$

(21)

Using Lemma 2, for the solution of the problem Eq. (20), we have:

$$ |v_{0} (x)|\,\, \le \,\beta^{ - 1} \mathop {\max }\limits_{0 \le s \le 1} |\phi (s)|\,\, $$

Thus, we obtain:

$$ |v_{0} (x)|\,\, \le \,\,C,\,\,\,\,\,0 \le x \le 1 $$

(22)

Applying the maximum principle to the problem Eq. (21), gives

$$ |v_{1} (x)|\,\, \le w(x)\,\,\, $$

(23)

where $w(x)$ is the solution of the following problem:

$$ \left\{ \begin{gathered} - \varepsilon w^{\prime\prime}(x) + \beta w(x) = 0,\,\,\,0 < x < 1, \hfill \\ w(0) = |v_{1} (0)|,\,\,\,w(1) = |v_{1} (1)| \hfill \\ \end{gathered} \right. $$

(24)

The solution of Eq. (4.24) has the form:

$$ w(x) \le C\left\{ {\frac{1}{\sqrt \varepsilon }\left( {e^{{ - \,\,\sqrt {\frac{\beta }{\varepsilon }} x}} + e^{{ - \sqrt {\frac{\beta }{\varepsilon }} (1 - x)}} } \right)} \right\}\, $$

(25)

Then, combining Eqs. (22), (23) and (25) we get the inequality,

$|u^{\prime}(x)|\,\, \le \,\,|v_{0} (x)| + |v_{1} (x)|\,\, \le C + w(x)\, = C + \frac{C}{\sqrt \varepsilon }\left( {e^{{ - \sqrt {\frac{\beta }{\varepsilon }} x}} + e^{{ - \sqrt {\frac{\beta }{\varepsilon }} (1 - x)}} } \right).$ Thus, the proof is completed. $\square$

Mesh generation

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

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

(26)

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

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

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

(27)

Formulation of the scheme

In order to formulate the scheme, consider the linear singularly perturbed differential equation in Eq. (12) subject to the conditions in Eq. (13) that gives the boundary value problem:

$$ \left\{ \begin{aligned}& - \varepsilon u^{\prime\prime}(x) + p(x)u(x) = f(x),\quad 0 < x < 1, \hfill \\ & u(0) = A, \hfill \\ & u(1) = \,B. \hfill \\ \end{aligned} \right. $$

(28)

Let us define the three the nodal points based finite difference approximation from the general multistep finite difference for the differential equation part as in the form:

$$ \sum\limits_{j = 0}^{2} {a_{j} U_{i - j + 1} = \sum\limits_{j = 0}^{2} {b_{j} U^{\prime\prime}_{i - j + 1} ,} } $$

(29)

where the coefficient parameters $a_{j}$ and $b_{j}$ are determined in terms the mesh parameter $h_{i}$. These parameters one can obtain in a similar way as on equidistant mesh. Hence, let the local truncation error estimated from Eq. (29) written as:

$$ T_{i} = a_{0} U_{i - 1} + a_{1} U_{i} + a_{2} U_{i + 1} - [b_{0} U^{\prime\prime}_{i - 1} + b_{1} U^{\prime\prime}_{i} + b_{2} U^{\prime\prime}_{i + 1} ]. $$

(30)

Assume that the function $u(x)$ has continuous derivatives of sufficiently fourth-order. Expanding the terms $U(x_{i \pm 1} )$ and $U^{\prime\prime}(x_{i \pm 1} )$ by Taylor’s series expansion about the point $x_{i}$ as:

$$ \left\{ \begin{aligned} & U_{i + 1} = U_{i} + h_{i + 1} U^{\prime}_{i} + \frac{{h_{i + 1}^{2} }}{2!}U^{\prime\prime}_{i} + \frac{{h_{i + 1}^{3} }}{3!}U^{\prime\prime\prime}_{i} + \frac{{h_{i + 1}^{4} }}{4!}U_{i}^{(4)} + \frac{{h_{i + 1}^{5} }}{5!}U_{i}^{(5)} + \frac{{h_{i + 1}^{6} }}{6!}U_{i}^{(6)} + O(h_{i + 1}^{7} ) \hfill \\ & U_{i - 1} = U_{i} - h_{i} U^{\prime}_{i} + \frac{{h_{i}^{2} }}{2!}U^{\prime\prime}_{i} - \frac{{h_{i}^{3} }}{3!}U^{\prime\prime\prime}_{i} + \frac{{h_{i}^{4} }}{4!}U_{i}^{(4)} - \frac{{h_{i}^{5} }}{5!}U_{i}^{(5)} + \frac{{h_{i}^{6} }}{6!}U_{i}^{(6)} + O(h_{i}^{7} ) \hfill \\ & U^{\prime\prime}_{i + 1} = U^{\prime\prime}_{i} + h_{i + 1} U^{\prime\prime\prime}_{i} + \frac{{h_{i + 1}^{2} }}{2!}U_{i}^{(4)} + \frac{{h_{i + 1}^{3} }}{3!}U_{i}^{(5)} + \frac{{h_{i + 1}^{4} }}{4!}U_{i}^{(6)} + O(h_{i + 1}^{7} ) \hfill \\ & U^{\prime\prime}_{i - 1} = U^{\prime\prime}_{i} - h_{i} U^{\prime\prime\prime}_{i} + \frac{{h_{i}^{2} }}{2!}U_{i}^{(4)} - \frac{{h_{i}^{3} }}{3!}U_{i}^{(5)} + \frac{{h_{i}^{4} }}{4!}U_{i}^{(6)} + O(h_{i}^{7} ) \hfill \\ \end{aligned} \right.\, $$

(31)

Then, substituting this Eq. (31) into Eq. (30) and grouping like terms gives:

$$ \begin{aligned} T_{i} & = (a_{0} + a_{1} + a_{2} )U_{i} + (a_{0} h_{i + 1} - a_{2} h_{i} )U^{\prime}_{i} \hfill \\ & \quad + (\frac{{a_{0} h_{i + 1}^{2} }}{2!} + \frac{{a_{2} h_{i}^{2} }}{2!} - (b_{0} + b_{1} + b_{2} ))U_{i}^{\prime \prime } \hfill \\ & \quad + (\frac{{a_{0} h_{i + 1}^{3} }}{3!} - \frac{{a_{2} h_{i}^{3} }}{3!} - (b_{0} h_{i + 1} - b_{2} h_{i} ))U^{\prime\prime\prime}_{i} \hfill \\ & \quad + (\frac{{a_{0} h_{i + 1}^{4} }}{4!} + \frac{{a_{2} h_{i}^{4} }}{4!} - (\frac{{b_{0} h_{i + 1}^{2} }}{2!} + \frac{{b_{2} h_{i}^{2} }}{2!}))U_{i}^{(4)} \hfill \\ & \quad + (\frac{{a_{0} h_{i + 1}^{5} }}{5!} - \frac{{a_{2} h_{i}^{5} }}{5!} - (\frac{{b_{0} h_{i + 1}^{3} }}{3!} - \frac{{b_{2} h_{i}^{3} }}{3!}))U_{i}^{(5)} \hfill \\ & \quad + (\frac{{a_{0} h_{i + 1}^{6} }}{6!} + \frac{{a_{2} h_{i}^{6} }}{6!} - (\frac{{b_{0} h_{i + 1}^{4} }}{4!} + \frac{{b_{2} h_{i}^{4} }}{4!}))U_{i}^{(6)} + ...\, \hfill \\ \end{aligned} $$

(32)

The method given in Eq. (29), is of order four if all the coefficients given in Eq. (32) are equal to zero except it is differ from zero after the coefficient of $U_{i}^{(6)} ,$ which gives the system of equation:

$$ \begin{aligned} & \left\{ \begin{aligned} & a_{0} + a_{1} + a_{2} = 0 \hfill \\ & a_{0} h_{i + 1} - a_{2} h_{i} = 0 \hfill \\ & a_{0} h_{i + 1}^{2} + a_{2} h_{i}^{2} = 2(b_{0} + b_{1} + b_{2} ) \hfill \\ & a_{0} h_{i + 1}^{3} - a_{2} h_{i}^{3} = 6(b_{0} h_{i + 1} - b_{2} h_{i} ) \hfill \\ & a_{0} h_{i + 1}^{4} + a_{2} h_{i}^{4} = 12(b_{0} h_{i + 1}^{2} + b_{2} h_{i}^{2} ) \hfill \\ & a_{0} h_{i + 1}^{5} - a_{2} h_{i}^{5} = 20(b_{0} h_{i + 1}^{3} - b_{2} h_{i}^{3} ) \hfill \\ & \end{aligned} \right.\, \hfill \\ & \quad a_{0} h_{i + 1}^{6} + a_{2} h_{i}^{6} \ne 30(b_{0} h_{i + 1}^{4} + b_{2} h_{i}^{4} ) \hfill \\ \end{aligned} $$

(33)

Adapting $b_{0} + b_{1} + b_{2} = 1,$ from the relation in uniform mesh, the solution of the system in Eq. (33) is determined by using the usually elimination method gives:

$$ \left\{ \begin{aligned} & a_{0} = \frac{2}{{h_{i + 1} (h_{i} + h_{i + 1} )}} \hfill \\ & a_{1} = \frac{ - 2}{{h_{i} h_{i + 1} }} \hfill \\ & a_{2} = \frac{2}{{h_{i} (h_{i} + h_{i + 1} )}} \hfill \\ \end{aligned} \right.\,\,\,\,\,\,\,\,\,\,and\,\,\,\,\,\,\,\,\,\,\, \left\{ \begin{aligned} & b_{0} = \frac{{h_{i + 1}^{2} + h_{i} h_{i + 1} - h_{i}^{2} }}{{6h_{i + 1} (h_{i} + h_{i + 1} )}} \hfill \\ & b_{1} = \frac{{4h_{i} h_{i + 1} (h_{i} + h_{i + 1} ) + h_{i}^{3} + h_{i + 1}^{3} }}{{6h_{i} h_{i + 1} (h_{i} + h_{i + 1} )}} \hfill \\ & b_{2} = \frac{{h_{i} h_{i + 1} + h_{i}^{2} - h_{i + 1}^{2} }}{{6h_{i} (h_{i} + h_{i + 1} )}} \hfill \\ \end{aligned} \right.\,\, $$

(34)

Using this we approximate the problem in Eq. (28) by Eq. (29), which can write in the form of:

$$ a_{0} U_{i - 1} + a_{1} U_{i} + a_{2} U_{i + 1} = b_{0} U^{\prime\prime}_{i - 1} + b_{1} U^{\prime\prime}_{i} + b_{2} U^{\prime\prime}_{i + 1} ,\,\, $$

(35)

Then, considering from the differential equation in Eq. (28), at the nodal point $x_{i} ,$ we have:

$$ U^{\prime\prime}_{i} = \frac{{p_{i} U_{i} - f_{i} }}{\varepsilon },\,\,\,\,\,U^{\prime\prime}_{i - 1} = \frac{{p_{i - 1} U_{i - 1} - f_{i - 1} }}{\varepsilon },\,\,\,\,\,\,{\text{and}}\,\,\,U^{\prime\prime}_{i + 1} = \frac{{p_{i + 1} U_{i + 1} - f_{i + 1} }}{\varepsilon }. $$

(36)

Substituting Eq. (36) into Eq. (35) and also using the values in Eq. (34), we obtain the three-term recurrence finite difference scheme

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

(37)

where

$$ \begin{aligned} E_{i} & = \frac{ - 2\varepsilon }{{h_{i + 1} (h_{i + 1} + h_{i} )}} + \frac{{p_{i - 1} (h_{i + 1}^{2} + h_{i} h_{i + 1} - h_{i}^{2} )}}{{6h_{i + 1} (h_{i + 1} + h_{i} )}},\quad F_{i} = \frac{2\varepsilon }{{h_{i} h_{i + 1} }} + \frac{{P_{i} (4h_{i} h_{i + 1} (h_{i + 1} + h_{i} ) + h_{i}^{3} + h_{i + 1}^{3} )}}{{6h_{i} h_{i + 1} (h_{i + 1} + h_{i} )}} \hfill \\ G_{i} & = \frac{ - 2\varepsilon }{{h_{i} (h_{i + 1} + h_{i} )}} + \frac{{p_{i + 1} (h_{i}^{2} + h_{i} h_{i + 1} - h_{i + 1}^{2} )}}{{6h_{i} (h_{i + 1} + h_{i} )}}, \hfill \\ H_{i} & = \left( {\frac{{h_{i + 1}^{2} + h_{i} h_{i + 1} - h_{i}^{2} }}{{6h_{i + 1} (h_{i + 1} + h_{i} )}}} \right)f_{i - 1} + \left( {\frac{{4h_{i} h_{i + 1} (h_{i + 1} + h_{i} ) + h_{i}^{3} + h_{i + 1}^{3} }}{{6h_{i} h_{i + 1} (h_{i + 1} + h_{i} )}}} \right)f_{i} + \left( {\frac{{h_{i}^{2} + h_{i} h_{i + 1} - h_{i + 1}^{2} }}{{6h_{i} (h_{i + 1} + h_{i} )}}} \right)f_{i + 1} . \hfill \\ \end{aligned} $$

Considering Eq. (4.2) and by solving this system of linear algebraic equations we obtain the approximate solution $U_{i} ,\,i = 0,1,2,...,N$ of $u(x_{i} )$ at the nodal points $x_{0} ,x_{1} ,x_{2} ,...,x_{N} .$

Convergence analysis

Let $u_{i}$ be the solution of Eq. (12) and $U_{i}$ be the solution to Eq. (37) at the nodal point $x_{i}$, then $z_{i} = u_{i} - U_{i}$,$0 \le i \le N$, with the estimate approximate error $z_{i}$, which satisfies the discrete problem

$$ \left\{ \begin{aligned} & Lz_{i} = R_{i} , \hfill \\ & z_{0} = 0 = z_{N,} \hfill \\ \end{aligned} \right. $$

(38)

where $R_{i}$ is the truncation error in Eq. (32).

Lemma 4

(Discrete Maximum Principle): Suppose that a mesh function $v_{i}$ satisfies $v_{0} \ge 0$ and $v_{N} \ge 0$. Then $L^{N} v_{i} \ge 0,$ $\forall \,\,\,1 \le i \le N - 1$ implies that $v_{i} \ge 0$, $\forall \,\,\,0 \le i \le N$.

Proof

Let $V_{i}$ be a value such that $V_{i}^{*} = \mathop {\min \,v_{i} }\limits_{1 \le i \le N - 1}$ and assume that $V_{i}^{*} < 0.$ Clearly, $i \in \left\{ {0,N} \right\}$ and therefore $V_{i}^{*\prime } = 0$ and $V_{i}^{*\prime \prime } \ge 0.$ Moreover, there is.

$$ \begin{gathered} L^{N} V_{i} = a_{0} V_{i + 1} + a_{1} V_{i} + a_{2} V_{i - 1} - [b_{0} V^{\prime\prime}_{i + 1} + b_{1} V^{\prime\prime}_{i} + b_{2} V^{\prime\prime}_{i - 1} ], \hfill \\ L^{N} V_{i} = a_{0} V_{i + 1} + a_{1} V_{i} + a_{2} V_{i - 1} - [b_{0} V^{\prime\prime}_{i + 1} + b_{1} V^{\prime\prime}_{i} + b_{2} V^{\prime\prime}_{i - 1} ] < 0. \hfill \\ \end{gathered} $$

which is a contradiction. It follows that $V_{i}^{*} \ge 0$ and thus $v_{i} \ge 0,\,\,\forall \,\,\,0 \le i \le N$.$\square$

Lemma 5

(Uniform Stability Estimate): If $U_{i}$ is any mesh function such that $U_{i} = 0 = U_{N}$, then

$$ |U_{i} |\, \le \beta^{ - 1} \mathop {\max }\limits_{1 \le i \le N - 1} |L^{N} U_{i} |,\,\,\,0 \le i \le N $$

(39)

Proof

Denote $Z_{i} = \beta^{ - 1} \mathop {\max }\limits_{1 \le i \le N - 1} |L^{N} U_{i} |,\,\,\,1 \le i \le N - 1.$ Introduce two mesh functions

$\xi_{i}^{ \pm } = Zp_{i} \pm U_{i}$, Clearly, $\xi_{0}^{ \pm } = 0 = \xi_{N}^{ \pm }$ and $\forall \,\,1 \le i \le N - 1$:

$$ L^{N} \xi_{i}^{ \pm } = Zp_{i} + L^{N} U_{i} \le 0. $$

Since $p_{i} \ge \beta \ge 0$ Lemma 4 implies that $\xi_{i}^{ \pm } \ge 0,\,\,\forall \,\,0 \le i \le N,$ and this completes the proof. $\square$

Lemma 6

(Error Boundness), The truncation error at the grid point $x_{i}$ is given by:

$$ L^{*} (u - U)_{i} = Lu_{i} - L^{N} U_{i} $$

(40)

From the formulated method, the non-zero estimated local truncation error provided in Eq. (32) with the conditions in Eqs. (33) and (34) written as:

$$ T_{i} = \left( {\frac{{a_{0} h_{i + 1}^{6} }}{6!} + \frac{{a_{2} h_{i}^{6} }}{6!} - \left( {\frac{{b_{0} h_{i + 1}^{4} }}{4!} - \frac{{b_{2} h_{i}^{4} }}{4!}} \right)} \right)U_{i}^{(6)} $$

(41)

where $a_{0} ,\,a_{2} ,\,b_{0} ,\,b_{2}$ are defined in Eq. (34).

Thus, from the relation $h_{i}^{4} > h_{i}^{6}$, we have

$$ |T_{i} |\, = \,|L^{*} (u - U)_{i} |\,\, \le \,\,C(h_{i}^{4} + h_{i + 1}^{4} ) $$

(42)

for $C = \frac{1}{24}\left\{ {||b_{0} ||_{\infty } + ||b_{2} ||_{\infty } } \right\}||U_{i}^{(6)} ||_{\infty }$ is arbitrary constant.

Let us consider $h = \mathop {\max }\limits_{\forall \,\,i \to \infty } \left\{ {h_{i} ,\,\,h_{i + 1} } \right\}$ and applying the uniform stability estimate (Lemma 5), yields:

$$ \mathop {\max }\limits_{0 \le i \le N} |u - U|\,\, \le Ch^{4} $$

(43)

Hence, the formulated method is fourth order convergent.

Numerical examples and results

To illustrate the applicability of the proposed method, we applied it model example. For the example whose exact solution is unknown, we use the double mesh principle to estimate the error and compute the experimental rate of convergence. The double mesh principle is given by:

$$ D_{\varepsilon }^{N} = \mathop {\max }\limits_{{x_{i} \in [0,1]^{N} }} |U_{i}^{N} - U_{i}^{2N} |,\,\, $$

(44)

where $U_{i}^{N}$ and $U_{i}^{2N}$ respectively, denotes the numerical solution obtained using N and 2N mesh intervals. Further, we calculate the order of convergence by the formula:

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

(45)

Example 1

Consider the semilinear singularly perturbed reaction–diffusion problem, [14]:

$$ \left\{ \begin{aligned} & - \varepsilon y^{\prime\prime}(x) - e^{{ - (x^{2} + y)}} = 0, \quad 0 < x < 1, \hfill \\ & y(0) = 0, \hfill \\ & y(1) = 1. \hfill \\ \end{aligned} \right. $$

The exact solution for this problem is unknown.

Example 2

Consider the linear singularly perturbed problem.

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

.

The source function is given by:

$\begin{aligned} f(x) &= 1 + x(1 - x) + (2\sqrt \varepsilon - x^{2} + x^{3} )\exp (\frac{ - (1 - x)}{{\sqrt \varepsilon }}) \hfill \\ & \quad + (2\sqrt \varepsilon - x(1 - x)^{2} )\exp (\frac{ - x}{{\sqrt \varepsilon }}) \hfill \\ \end{aligned}$.

The exact solution is $y(x) = 1 + (x - 1)\exp (\frac{ - x}{{\sqrt \varepsilon }}) - x\exp (\frac{ - (1 - x)}{{\sqrt \varepsilon }}).$

Conclusion

In this paper, fourth-order fitted mesh scheme is presented for solving semilinear singularly perturbed reaction–diffusion problem. From this problem, the nonlinear part is linearized by the quasilinearization technique. The convergence analysis of the described method established theoretically as well as confirmed in numerical illustration that is fourth order convergent. To validate the applicability of method, model examples are considered and numerical results investigated in tabular and graphic forms. Specifically, results are expressed in terms of maximum absolute errors, and rate of convergence. The result in Table 1 indicates that the comparison of maximum absolute errors for the proposed method and the methods in [14]. Further, we observe that as the number mesh increases the maximum absolute error decreases in each row implying that the proposed method is convergent. Tables 2 and 3 indicates that the rate of convergence for the described method is fourth-order convergent. This confirms the theoretical investigations given by Eq. (43). Also, Fig. 1 used to visualize the effect of perturbation parameter and the boundary layer behaviors. Furthermore, the main originality of the suggested method describe and examined in terms the obtained more accurate solutions with higher order of convergence as our original contributions. Generally, the method is fourth-order convergent and gives more accurate solution than some existing methods in the literature.

Table 1 Comparison of maximum absolute errors for Example 1

Table 2 Computed rate of convergence for Example 1

Table 3 Computed maximum absolute errors and rate of convergence for Example 2

Fig. 1

Numerical solution for Example 1 with $N = 64$, $\varepsilon = 10^{ - 2}$ and $\varepsilon = 10^{ - 5}$ respectively

Limitations

During the quasilinearization process, it was fixed to the first iteration. If more number of iterations were done, then the scheme can have more accurate solution than the presented results. Additional, the scheme can more illustrate the physical behaviour of the problem under consideration.