Fourth-order stable central difference with Richardson extrapolation method for second-order self-adjoint singularly perturbed boundary value problems
Muslima Kedir Siraj1,
Gemechis File Duressa1 &
Tesfaye Aga Bullo1
This study introduces a stable central difference method for solving second-order self-adjoint singularly perturbed boundary value problems. First, the solution domain is discretized. Then, the derivatives in the given boundary value problem are replaced by finite difference approximations and the numerical scheme that provides algebraic systems of equations is developed. The obtained system of algebraic equations is solved by Thomas algorithm. The consistency and stability that guarantee the convergence of the scheme are investigated. The established convergence of the scheme is further accelerated by applying the Richardson extrapolation which yields sixth order convergent. To validate the applicability of the method, two model examples are solved for different values of perturbation parameter ε and different mesh size h. The proposed method approximates the exact solution very well. Moreover, the present method is convergent and gives more accurate results than some existing numerical methods reported in the literature.
Any differential equation obtained from a given differential equation and having the property that its solution is an integrating factor of the other is known as adjoint differential equation. Self-adjoint singularly perturbed differential equation is a differential equation whose highest order derivative is multiplied by a small positive parameter and that has the same solution as its adjoint equation [1, 2]. In a singularly perturbed problem, small positive parameter affects the highest order derivative(s) of the differential equation which gives rise to large gradients in the solution over narrow regions of the domain, so that the presence of a small perturbation parameter in the differential equation typically leads to boundary layers in the solution, which makes the convergence analysis very difficult [3]. As Miller et al. [2], boundary layer is a region of the independent variable over which the dependent variable changes rapidly.
Singularly perturbed second-order boundary value problem occur very frequently in fluid motion, chemical reactor theory, elasticity, diffusion in polymer, reaction-diffusion equation, control of chaotic system, and so on [4]. Upon setting ε = 0, if the order of singularly perturbed differential equations is reduced by one, then the problem is called convection-diffusion type and if the order is reduced by two, it is called reaction-diffusion type. Hence, second-order singularly perturbed self-adjoint ordinary differential equations are types of reaction-diffusion problems. Since the solution of this problem exhibits one or two layers, the existing numerical methods give good results only when the mesh size h is smaller than the perturbation parameter ε (i.e., h ≤ ε). But it is an expensive and time-consuming process. If we take h ≥ ε, the existing classical numerical methods produce oscillatory solution and pollute the solution in the entire interval, because of boundary layer behavior. In connection to this, there are some numerical methods suggested by various authors for solving self-adjoint singular perturbation problems, namely, initial value technique [5], quintic spline method [6], non-polynomial spline functions method [7], difference scheme using cubic spline [8], finite difference method with variable mesh [9], fitted mesh B-spline collocation method [10], higher order numerical methods [11].
More recently, Fasika et al. [12,13,14] and Feyisa and Gemechis [15] have developed a higher (fourth, sixth, eighth, and tenth) order compact finite difference method to solve singularly perturbed reaction-diffusion problems. These authors developed higher order compact finite difference methods for the constant coefficients of diffusion and reaction terms of the problem. Even though their methods give more accurate numerical solutions, it is restricted to treat constant coefficient problem. Also, other scholars, Terefe et al. [16] and Yitbarek et al. [17], have presented fourth- and sixth-order stable central difference method, respectively, for solving self-adjoint singularly perturbed two-point boundary value problem. Therefore, the main objective of this study is to develop a stable and more accurate numerical method that works for solving both constant and variable coefficient second-order self-adjoint singularly perturbed boundary value problems.
In this paper, we planned a fourth-order stable central difference with Richardson extrapolation method for solving second-order self-adjoint singularly perturbed boundary value problems. First, the derivative in the given differential equation is replaced by the finite difference approximations. Then, the ordinary differential equation converts to a linear system of algebraic equations, and these algebraic equations are transformed to a tri-diagonal system, which can easily be solved by the Thomas algorithm. Further, coding of the program in MATLAB software for the obtained tri-diagonal system has been performed. To validate the applicability of the method, some model examples are considered for numerical experimentation. Both the theoretical and numerical rates of convergence of the scheme have been investigated.
Formulation of the method
Consider the singularly perturbed self-adjoint boundary value problem of the form:
$$ -\varepsilon {\left(a(x){y}^{\prime }(x)\right)}^{\prime }+b(x)y(x)=g(x),\kern0.96em x\in \varOmega := \left(0,1\right) $$
subject to the boundary conditions:
$$ y(0)=\alpha $$
$$ y(1)=\beta $$
where ε is a perturbation parameter that satisfies 0 < ε < < 1, α, β are given constants and the functions a(x) ≠ 0, b(x) ≠ 0 and g(x) are assumed to be sufficiently continuous differentiable functions. By product rule differentiation, Eq. (1) can be re-written as:
$$ -\varepsilon {y}^{{\prime\prime} }(x)+p(x){y}^{\prime }(x)+q(x)y(x)=f(x) $$
$$ p(x)=\frac{-\varepsilon {a}^{\prime }(x)}{a(x)},\kern1.2em q(x)=\frac{b(x)}{a(x)} $$
$$ f(x)=\frac{g(x)}{a(x)}. $$
In order to develop the finite difference method, the interval [0, 1] is divided into N equal sub-intervals with set of grid points xi = x0 + ih, for i = 0, 1, 2, ..., N, where \( h=\frac{1}{N} \). For convenience, let \( p\left({x}_i\right)={p}_i,q\left({x}_i\right)={q}_i,y\left({x}_i\right)={y}_i,{y}^{\prime}\left({x}_i\right)={y}_i^{\prime },...,{y}^{(n)}\left({x}_i\right)={y}_i^{(n)}. \)
Assume that y(x) has continuous higher order derivatives on [0, 1], and to develop the fourth-order stable central difference scheme, we use Taylor's series expansion in order to get central difference formula for yi′′ andyi′.
$$ {y}_{i+1}={y}_i+h{y_i}^{\prime }+\frac{h^2}{2!}{y_i}^{\prime \prime }+\frac{h^3}{3!}{y_i}^{\prime \prime \prime }+\frac{h^4}{4!}{y}_i^{(4)}+\frac{h^5}{5!}{y}_i^{(5)}+\frac{h^6}{6!}{y}_i^{(6)}+... $$
$$ {y}_{i-1}={y}_i-h{y_i}^{\prime }+\frac{h^2}{2!}{y_i}^{\prime \prime }-\frac{h^3}{3!}{y_i}^{\prime \prime \prime }+\frac{h^4}{4!}{y}_i^{(4)}-\frac{h^5}{5!}{y}_i^{(5)}+\frac{h^6}{6!}{y}_i^{(6)}+... $$
From Eqs. (4) and (5), we have:
$$ {y_i}^{\prime }=\frac{y_{i+1}-{y}_{i-1}}{2h}-\frac{h^2}{6}{y}_i^{{\prime\prime\prime} }-\frac{h^4}{120}{y}_i^{(5)}+{\tau}_1\kern0.48em \mathrm{and}\kern0.36em {y_i}^{\prime \prime }=\frac{y_{i+1}-2{y}_i+{y}_{i-1}}{h^2}-\frac{h^2}{12}{y}_i^{(4)}+{\tau}_2 $$
where \( {\tau}_1=\frac{-{h}^6}{7!}{y_i}^{(7)} \) and \( {\tau}_2=\frac{-{h}^4}{360}{y_i}^{(6)} \)
Substituting Eq. (6) into the discrete form of Eq. (3) gives:
$$ {q}_i{y}_i+\frac{p_i}{2h}\left({y}_{i+1}-{y}_{i-1}\right)-\frac{\varepsilon }{h^2}\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)-\frac{p_i{h}^2}{6}{y}_i^{{\prime\prime\prime} }+\frac{\varepsilon {h}^2}{12}{y}_i^{(4)}-\frac{p_i{h}^4}{120}{y}_i^{(5)}+{\tau}_0={f}_i $$
$$ {\tau}_0={p}_i{\tau}_1-\varepsilon {\tau}_2 $$
Differentiating Eq. (3) successively and considering at the nodal points yields:
$$ {y_i}^{\prime \prime \prime }=\frac{1}{\varepsilon}\left({p}_i{y_i}^{\prime \prime }+\left({p_i}^{\prime }+{q}_i\right){y_i}^{\prime }+{q_i}^{\prime }{y}_i-{f_i}^{\prime}\right) $$
$$ {y_i}^{(4)}=\frac{1}{\varepsilon}\left({p}_i{y_i}^{\prime \prime \prime }+\left(2{p_i}^{\prime }+{q}_i\right){y_i}^{\prime \prime }+\left({p_i}^{\prime \prime }+2{q_i}^{\prime}\right){y}_i^{\prime }+{q_i}^{\prime \prime }{y}_i-{f_i}^{\prime \prime}\right) $$
$$ {y_i}^{(5)}=\frac{1}{\varepsilon}\left({p}_i{y_i}^{(4)}+\left(3{p}_i^{\prime }+{q}_i\right){y}_i^{{\prime\prime\prime} }+\left(3{p_i}^{\prime \prime }+3{q_i}^{\prime}\right){y_i}^{\prime \prime }+\left({p}_i^{{\prime\prime\prime} }+3{q}^{{\prime\prime}}\right){y}_i^{\prime }+q{i}^{{\prime\prime\prime} }{y}_i-{f}_i^{{\prime\prime\prime}}\right) $$
Using Eq. (10), the term which contains yi(5) from Eq. (7) becomes:
$$ {\displaystyle \begin{array}{l}\frac{-{p}_i{h}^4}{120}{y_i}^{(5)}=-\frac{{p_i}^2{h}^4}{120\varepsilon }{y_i}^{(4)}-\frac{p_i{h}^{(4)}}{120\varepsilon}\left(3{p_i}^{\prime }+{q}_i\right){y_i}^{\prime \prime \prime }-\frac{p_i{h}^{(4)}}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}_i^{\prime}\right){y}_i^{{\prime\prime}}\\ {}\kern3.719998em \frac{-{p}_i{h}^4}{120\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right){y}_i^{\prime }-\frac{p_i{q_i}^{\prime \prime \prime }{h}^4}{120\varepsilon }{y}_i+\frac{p_i{h}^4}{120\varepsilon }{f}_i^{{\prime\prime\prime}}\end{array}} $$
Also, from Eqs. (4) and (5), we have the central finite difference approximation:
$$ {y_i}^{\prime }=\frac{y_{i+1}-{y}_{i-1}}{2h}+{\tau}_3 $$
$$ {y_i}^{\prime \prime }=\frac{y_{i+1}-2{y}_i+{y}_{i-1}}{h^2}+{\tau}_4 $$
$$ {\tau}_3=\frac{-{h}^2}{6}{y}_i^{{\prime\prime\prime} } $$
$$ {\tau}_4=-\frac{h^2}{12}{y}_i^{(4)} $$
Putting Eq. ((12), into Eq. (11) gives:
$$ {\displaystyle \begin{array}{l}\frac{-{p}_i{h}^4}{120}{y_i}^{(5)}=-\frac{{p_i}^2{h}^4}{120\varepsilon }{y_i}^{(4)}-\frac{p_i{h}^{(4)}}{120\varepsilon}\left(3{p_i}^{\prime }+{q}_i\right){y_i}^{\prime \prime \prime}\\ {}\kern3.479999em -\left(\frac{p_i{h}^2}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}_i^{\prime}\right)\right)\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)\\ {}\kern3.359999em -\left(\frac{p_i{h}^3}{120\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right)\right)\left({y}_{i+1}-{y}_{i-1}\right)-\frac{p_i{q_i}^{\prime \prime \prime }{h}^4}{120\varepsilon }{y}_i+\frac{p_i{h}^4}{120\varepsilon }{f}_i^{{\prime\prime\prime} }+{\tau}_5\end{array}} $$
$$ {\tau}_5=-\frac{p_i{h}^4}{120\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right){\tau}_3-\frac{p_i{h}^4}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}^{\prime}\right){\tau}_4. $$
Substituting Eq. (13) into Eq. (7) and rearranging, we get:
$$ {\displaystyle \begin{array}{l}\left({q}_i-\frac{p_i{q}_i^{{\prime\prime\prime} }{h}^4}{120\varepsilon}\right){y}_i+\left(\frac{p_i}{2h}-\frac{p_i{h}^3}{240\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right)\right)\left({y}_{i+1}-{y}_{i-1}\right)\\ {}-\left(\frac{\varepsilon }{h^2}+\frac{p_i^2}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}_i^{\prime}\right)\right)\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)-\left(\frac{pi{h}^2}{6}+\frac{pi{h}^4}{120\varepsilon}\left(3{p}_i^{\prime }+ qi\right)\right){y}^{{\prime\prime\prime}}\\ {}+\left(\frac{\varepsilon {h}^2}{12}-\frac{p_i^2{h}^2}{120\varepsilon}\right){y}_i^4+{\tau}_6={f}_i-\frac{p_i{h}^4}{120\varepsilon }{f}_i^{{\prime\prime\prime}}\end{array}} $$
$$ {\tau}_6={\tau}_0+{\tau}_5 $$
Again, using Eq. (9), the term which contains \( {y}_i^{(4)} \) from Eq. (14) becomes:
$$ {\displaystyle \begin{array}{l}\left(\frac{\varepsilon {h}^2}{12}-\frac{p_i^2{h}^4}{120\varepsilon}\right){y_i}^4=\left(\frac{p_i{h}^2}{12}-\frac{p_i^3{h}^4}{120{\varepsilon}^2}\right){y}_i^{{\prime\prime\prime} }+\left(\frac{q_i^{{\prime\prime} }{h}^2}{12}-\frac{p_i^2{q}_i^{{\prime\prime} }{h}^4}{120{\varepsilon}^2}\right){y}_i+{\tau}_7\\ {}\kern4.919997em +\left(\left(\frac{1}{12}-\frac{p_i^2{h}^2}{120{\varepsilon}^2}\right)\left(2{p}_i^{\prime }+ qi\right)\right)\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)\\ {}\kern4.199998em +\left(\left(\frac{h}{24}-\frac{p_i^2{h}^3}{240{\varepsilon}^2}\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right)\right)\left({y}_{i+1}-{y}_{i-1}\right)-\left(\frac{h^2}{12}-\frac{p_i^2{h}^4}{120{\varepsilon}^2}\right){f}_i^{{\prime\prime}}\end{array}} $$
where \( {\tau}_7=\left(\frac{h^2}{12}-\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\right)\left(2{p}_i^{\prime }+{q}_i\right){\tau}_4+\left(\frac{h^2}{12}-\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right){\tau}_3 \)
Substituting Eq. (15) into Eq. (14) yields:
$$ {\displaystyle \begin{array}{l}\left({q}_i-\frac{p_i{q}_i^{{\prime\prime\prime}}\;{h}^4}{120\varepsilon }+\frac{q_i^{{\prime\prime\prime}}\;{h}^2}{12}-\frac{p_i^2{q}_i^{{\prime\prime\prime}}\;{h}^4}{120{\varepsilon}^2}\right){y}_i\\ {}+\left(\frac{p_i}{2h}-\frac{p_i{h}^3}{240\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right)+\left(\frac{h}{24}-\frac{p_i^2{h}^3}{240{\varepsilon}^2}\right)\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right)\left({y}_{i+1}-{y}_{i-1}\right)\\ {}-\left(\frac{\varepsilon }{h^2}+\frac{p_i{h}^2}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}_i^{\prime}\right)-\left(\frac{1}{12}-\frac{p_i^2{h}^2}{120{\varepsilon}^2}\right)\right)\left(2{p}_i^{\prime }+{q}_i\right)\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)\\ {}+\left(\frac{p_i{h}^2}{12}+\frac{p_i{h}^4}{120\varepsilon}\left(3{p}_i^{\prime }+ qi\right)-\frac{p_i{h}^2}{6}+\frac{{p_i}^3{h}^4}{120{\varepsilon}^2}\;\right){y}^{{\prime\prime\prime} }+{\tau}_8\\ {}\kern7.799994em ={f}_i+\left(\frac{h^2}{12}+\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\;\right){f}_i^{{\prime\prime} }-\frac{p_i{h}^4}{120\varepsilon }{f}_i^{{\prime\prime\prime}}\end{array}} $$
$$ {\tau}_8={\tau}_6+{\tau}_7. $$
For simplicity, let
\( {A}_i={q}_i-\frac{p_i{q}_i^{{\prime\prime\prime} }{h}^4}{120\varepsilon }+\frac{q_i^{{\prime\prime\prime} }{h}^2}{12}-\frac{p_i^2{q}_i^{{\prime\prime\prime} }{h}^4}{120{\varepsilon}^2}, \) \( {B}_i=\frac{p_i}{2h}-\frac{p_i{h}^3}{240\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right)+\left(\frac{h}{24}-\frac{p_i^2{h}^3}{240{\varepsilon}^2}\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right) \)
$$ {C}_i=\frac{\varepsilon }{h^2}+\frac{p_i{h}^2}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}_i^{\prime}\right)-\left(\frac{1}{12}-\frac{p_i^2{h}^2}{120{\varepsilon}^2}\right)\left(2{p}_i^{\prime }+{q}_i\right) $$
\( {D}_i=\frac{p_i{h}^2}{12}-\frac{p_i{h}^4}{120\varepsilon}\left(3{p}_i^{\prime }+{q}_i\right)-\frac{p_i{h}^2}{6}-\frac{{p_i}^3{h}^4}{120{\varepsilon}^2}, \) \( Hh(i)={f}_i+\left(\frac{h^2}{12}+\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\;\right){f}_i^{{\prime\prime} }-\frac{p_i{h}^4}{120\varepsilon }{f}^{{\prime\prime\prime} } \)
Then, Eq. (16) re-written as:
$$ {A}_i{y}_i+{B}_i\left({y}_{i+1}-{y}_{i-1}\right)-{C}_i\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)+{D}_i{y}_i^{{\prime\prime\prime} }+{\tau}_8= Hh(i) $$
Lastly, using Eq. (8), the term that contains \( {y}_i^{{\prime\prime\prime} } \) from Eq. (17) becomes:
$$ {D}_i{y_i}^{\prime \prime \prime }=\frac{p_i{D}_i}{\varepsilon}\left({y}_{i+1}-2{y}_i+{y}_{i-1}\right)+\frac{D_i}{2 h\varepsilon}\left({p_i}^{\prime }+{q}_i\right)\left({y}_{i+1}-{y}_{i-1}\right)+\frac{D_i{q_i}^{\prime }}{\varepsilon }{y}_i-\frac{D_i}{\varepsilon }{f_i}^{\prime }+{\tau}_9 $$
$$ {\tau}_9=\frac{p_i{D}_i}{\varepsilon }{\tau}_4+\frac{D_i}{\varepsilon }{\tau}_3 $$
Putting Eq. (18) into Eq. (17), and write in three-term recurrence relation:
$$ -{E}_i{y}_{i-1}+{F}_i{y}_i-{G}_i{y}_{i+1}+{\tau}_{10}={H}_i $$
where\( {E}_i={C}_i-\frac{p_i{D}_i}{\varepsilon {h}^2}+{B}_i+\frac{D_i}{2 h\varepsilon}\left({p_i}^{\prime }+{q}_i\right), \) \( {F}_i={A}_i+\frac{D_i{q}_i^{\prime }}{\varepsilon }+2\left({C}_i-\frac{P_i{D}_i}{\varepsilon {h}^2}\right) \)\( {G}_i={C}_i-\frac{P_i{D}_i}{\varepsilon {h}^2}-{B}_i-\frac{D_i}{2 h\varepsilon}\left({p_i}^{\prime }+{q}_i\right) \) and \( {H}_i= Hh(i)+\frac{D_i}{\varepsilon }{f}_i^{\prime } \)with truncation error: τ10 = τ8 + τ9.
Richardson extrapolation
The basic idea behind extrapolation is that whenever the leading term in the error for an approximation formula is known, we can combine two approximations obtained from the formula using different values of the mesh sizes h and 0.5h to obtain a higher order approximation and the technique is known as Richardson extrapolation. This procedure is a convergence acceleration technique which consists of a linear combination of two computed approximations of a solution (applied on two nested meshes). The linear combination turns out to be a better approximation.
Since the truncation error of the formulated method Eq. (19) is O(h4), we have
$$ \mid y\left({x}_i\right)-{Y}_N\mid \le C\left({h}^4\right) $$
where y(xi) and YN are exact and approximate solutions respectively, C is constant independent of mesh sizes h.
Let Ω2N be the mesh obtained by bisecting each mesh interval in ΩN and denote the approximation of the solution on Ω2N by Y2N. Consider Eq. (20) works for any h ≠ 0, which implies:
$$ y\left({x}_i\right)-{Y}_N\le C\left({h}^4\right)+{R}^N,\kern0.96em {x}_i\in {\varOmega}^N $$
So that it also works for any \( \frac{h}{2}\ne 0 \) and yields:
$$ y\left({x}_i\right)-{Y}_{2N}\le C\left({\left(\frac{h}{2}\right)}^4\right)+{R}^{2N},\kern0.96em {x}_i\in {\varOmega}^{2N} $$
where the remainders, RN and R2N, are of O(h6). A combination of inequalities in Eqs. (21) and (22) leads to 15y(xi) − (16Y2N − YN) ≈ O(h6), which suggests that
$$ {\left({Y}_N\right)}^{ext}=\frac{1}{15}\left(16{Y}_{2N}-{Y}_N\right) $$
is also an approximation of y(xi). Using this approximation to evaluate the truncation error, we obtain:
$$ \mid y\left({x}_i\right)-{\left({Y}_N\right)}^{ext}\mid \le C{h}^6 $$
Now, using the solutions obtained by the scheme given by Eq. (19), we get another third solution in terms of the two by Eq. (23). This is Richardson extrapolation method for the fourth-order finite-difference scheme only to accelerate the rate of convergence to sixth order.
Consistency of the method
Local truncation errors refer to the differences between the original differential equation and its finite difference approximations at grid points. Local truncation errors measure how well a finite difference discretization approximates the differential equation [18]. In our case, the last truncation error in Eq. (19) is τ10 = τ8 + τ9.
But, from Eqs. (16) and (17), we have τ8 = τ6 + τ7and \( {\tau}_9=\frac{p_i{D}_i}{\varepsilon }{\tau}_4+\frac{D_i}{\varepsilon }{\tau}_3 \), So that:
$$ {\tau}_{10}={\tau}_6+{\tau}_7+\frac{p_i{D}_i}{\varepsilon }{\tau}_4+\frac{D_i}{\varepsilon }{\tau}_3 $$
$$ {\tau}_{10}={\tau}_0+{\tau}_5+{\tau}_7+\frac{p_i{D}_i}{\varepsilon }{\tau}_4+\frac{D_i}{\varepsilon }{\tau}_3, $$
because of Eq. (13)
Also, from Eqs. (13) and (15), we have:
$$ {\displaystyle \begin{array}{l}{\tau}_5=-\frac{p_i{h}^4}{120\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right){\tau}_3-\frac{p_i{h}^4}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}^{\prime}\right){\tau}_4\\ {}{\tau}_7=\left(\frac{h^2}{12}-\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\right)\left(2{p}_i^{\prime }+{q}_i\right){\tau}_4+\left(\frac{h^2}{12}-\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right){\tau}_3\end{array}} $$
Hence, the truncation errors are re-written as:
$$ {\displaystyle \begin{array}{l}{\tau}_{10}={\tau}_0+\left(\frac{D_i}{\varepsilon }+\left(\frac{h^2}{12}-\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right)-\frac{p_i{h}^4}{120\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right)\right){\tau}_3\\ {}\kern2.04em +\left(\frac{p_i{D}_i}{\varepsilon }+\left(\frac{h^2}{12}-\frac{{p_i}^2{h}^4}{120{\varepsilon}^2}\right)\left(2{p}_i^{\prime }+{q}_i\right)-\frac{p_i{h}^4}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}^{\prime}\right)\right){\tau}_4\end{array}} $$
Again, from Eqs. (6), (7), and ((12) into Eq. (25) and after rearranging yields:
$$ \mid TE\mid \le C{h}^4 $$
where TE = τ10 and
$$ {\displaystyle \begin{array}{l}C=\mid \frac{1}{6}{y}_i^{{\prime\prime\prime}}\;\left(\frac{p_i}{12\varepsilon }+\frac{p_i{h}^2}{120{\varepsilon}^2}\left(3{p}_i^{\prime }+ qi\right)+\frac{{p_i}^3{h}^2}{120{\varepsilon}^3}+\left(\frac{1}{12}-\frac{{p_i}^2{h}^2}{120{\varepsilon}^2}\right)\left({p}_i^{{\prime\prime} }+2{q}_i^{\prime}\right)-\frac{p_i{h}^2}{120\varepsilon}\left({p}_i^{{\prime\prime\prime} }+3{q}_i^{{\prime\prime}}\right)\right)\\ {}\kern0.72em +\frac{1}{12}{y}_i^{(4)}\left(\frac{p_i^2}{12\varepsilon }+\frac{p_i^2{h}^2}{120{\varepsilon}^2}\left(3{p}_i^{\prime }+ qi\right)+\frac{{p_i}^4{h}^2}{120{\varepsilon}^3}+\left(\frac{1}{12}-\frac{{p_i}^2{h}^2}{120{\varepsilon}^2}\right)\left(2{p}_i^{\prime }+{q}_i\right)-\frac{p_i{h}^2}{120\varepsilon}\left(3{p}_i^{{\prime\prime} }+3{q}^{\prime}\right)\right)\\ {}\kern1.32em +\frac{\varepsilon }{360}{y_i}^{(6)}-\frac{p_i{h}^2}{5040}{y_i}^{(7)}\mid \\ {}\kern1.56em \end{array}} $$
Thus, the developed scheme without applying Richardson extrapolation is fourth order accurate or order of convergence is O(h4). As Zhilin et al. [18], a finite difference scheme is called consistent if the limit of truncation error (TE) is equal to zero as the mesh size h goes to zero. Hence, this definition of consistency on the proposed method which is given in Eq. (19) with the local truncation error in Eqs. (24) and (26) satisfied as:
$$ \underset{h\to 0}{\lim } TE=\underset{h\to 0}{\lim }C{h}^4=\underset{h\to 0}{\lim }C{h}^6=0 $$
Thus, the proposed method is consistent.
Stability of the method
Consider the developed scheme in Eq. (19) which is given by:
$$ -{E}_i{y}_{i-1}+{F}_i{y}_i-{G}_i{y}_{i+1}={H}_i $$
But, the coefficients Ei, Fi and Gi are given in terms of Ai, Bi, Ci and Di with its values stated in Eq. (17). If we multiply both sides of Eq. (17) by h2 and consider the limit ash → 0, we get:
$$ {A}_i={B}_i={D}_i=0\;\mathrm{and}\;{C}_i=\varepsilon $$
Using the values in Eq. (27), the coefficients Ei, Fi, and Giin Eq. (19) becomes:
$$ {E}_i={G}_i=\varepsilon \kern0.48em \mathrm{and}\kern0.24em {F}_i=2\varepsilon $$
Considering both Eqs. (27) and (28), into Eq. (19), which can be written in matrix form:
$$ MY=H $$
where the matrices:
\( M=\left[\begin{array}{cccccc}2\varepsilon & -\varepsilon & 0& \cdots & \cdots & 0\\ {}-\varepsilon & 2\varepsilon & -\varepsilon & 0& \vdots & \vdots \\ {}0& -\varepsilon & 2\varepsilon & \ddots & 0& 0\\ {}0& 0& \ddots & \ddots & -\varepsilon & 0\\ {}\vdots & \vdots & 0& -\varepsilon & 2\varepsilon & -\varepsilon \\ {}0& \cdots & 0& 0& -\varepsilon & 2\varepsilon \end{array}\right] \), \( Y=\left[\begin{array}{c}{y}_1\\ {}{y}_2\\ {}\vdots \\ {}\vdots \\ {}{y}_{N-2}\\ {}{y}_{N-1}\end{array}\right]\kern0.36em \mathrm{and}\kern0.6em H=\left[\begin{array}{c}{h}^2{H}_1+\varepsilon {y}_0\\ {}{h}^2{H}_2\\ {}\vdots \\ {}\vdots \\ {}{h}^2{H}_{N-2}\\ {}{h}^2{H}_{N-1}+\varepsilon {y}_N\end{array}\right] \)
Here, M is a tri-diagonal matrix. M is irreducible if its co-diagonals contain non-zero elements only. The co-diagonal contains Ei, Gi. It is easily seen that, for sufficiently small h(i. e. h → 0), Ei ≠ 0 and Gi ≠ 0, ∀i = 1, 2, ⋯, N − 1.
Hence, M is irreducible. Again, one can observe that ∣Ei ∣ > 0 and ∣ Gi ∣ > 0 and in each row of M, the sum of the two off-diagonal elements is less than or equal to the modulus of the diagonal element (i.e., ∣Fi ∣ ≥ ∣ Ei ∣ + ∣ Gi∣). This proves the diagonal dominance of M. Under these conditions, the Thomas algorithm is stable for sufficiently small h, as shown in [19].
As proved by Smith [20], the eigenvalues of a tri-diagonal matrix (N − 1) × (N − 1) of matrix M are:
$$ {\lambda}_s={F}_i-2\sqrt{E_i{G}_i}\cos \frac{s\pi}{N}=2\varepsilon \left(1-\cos \frac{s\pi}{N}\right),\kern0.72em s=1,2,\dots, N-1 $$
Also, from trigonometric identity, we have \( 1-\cos \frac{s\pi}{N}=2{\sin}^2\frac{s\pi}{2N} \). Hence, the eigenvalues of matrix M can be re-written as:
$$ {\lambda}_s=2\varepsilon \left(2{\sin}^2\frac{s\pi}{2N}\right)=4\varepsilon {\sin}^2\frac{s\pi}{2N}\kern0.48em \le 4\varepsilon $$
A finite difference method for the BVPs is stable if M is invertible and
$$ \left\Vert {M}^{-1}\right\Vert \le C,\kern1em \forall 0<h<{h}_0 $$
where C and h0are two constants that are independent of h,[20].
Since, matrix M is symmetric also its inverse matrix M−1 is symmetric and the eigenvalues M−1 is given by \( \frac{1}{\lambda_s} \), we have
\( \left\Vert {M}^{-1}\right\Vert =\frac{1}{\lambda_s}=\frac{1}{4\varepsilon}\le C, \) where C is independent of h.
Thus, the developed scheme in Eq. (19) is stable.
A consistent and stable finite difference method is convergent by Lax's equivalence theorem [20]. Hence, as we have shown above, the proposed method is satisfying the criteria for both consistency and stability which are equivalents to convergence of the method.
Numerical examples and results
In order to test the validity of the proposed method and to demonstrate their convergence computationally, we have taken two model examples of singularly perturbed self-adjoint second-order two-point boundary value problems with exact solutions. The maximum absolute errors (AE) at the nodal points are given by \( \mid AE\mid =\underset{1\le i\le N-1}{\max}\mid y\left({x}_i\right)-{\left({Y}_N\right)}^{ext}\mid \). And the rate of convergence (R) can be calculated by the formula:
$$ R=\frac{\log {\left({Y}_N\right)}^{ext}-\log {\left({Y}_{2N}\right)}^{ext}}{\log 2} $$
where y(xi) and (YN)ext are exact solution and numerical solution, respectively, at the nodal point xi. And for the rate of convergence YN and Y2N are the numerical solutions obtained on the mesh size h and \( \frac{h}{2} \), respectively.
Consider the singularly perturbed self-adjoint problem:
$$ -\varepsilon {\left(\left(1+{x}^2\right){y}^{\prime }(x)\right)}^{\prime }+\left(1+x-{x}^2\right)y(x)=f(x),\kern0.72em 0<x<1 $$
subject to the boundary conditions y(0) = y(1) = 0, where f (x) is chosen such that the exact solution is given by: \( y(x)=1+\left(x-1\right)\;{e}^{\frac{-x}{\sqrt{\varepsilon }}}-x{e}^{\frac{\left(1-x\right)}{\sqrt{\varepsilon }}} \)
Consider the following self-adjoint singular perturbation problem:
$$ -\varepsilon {y}^{{\prime\prime} }(x)+\frac{4}{{\left(x+1\right)}^4}\left(\left(x+1\right)\sqrt{\varepsilon}\right)y(x)=f(x),\kern0.96em 0<x<1 $$
with boundary conditions y(0) = 2 and y(1) = − 1, where f(x) is chosen such that the exact solution is given by: \( y(x)=-\cos \left(\frac{4\pi x}{x+1}\right)+\frac{3\left({\exp}^{\left(\frac{-2x}{\sqrt{\varepsilon}\left(x+1\right)}\right)}-{\exp}^{\left(\frac{-1}{\sqrt{\varepsilon }}\right)}\right)}{1-{\exp}^{\left(\frac{-1}{\sqrt{\varepsilon }}\right)}} \)
In this paper, we described the fourth-order stable central difference method with Richardson extrapolation for solving second-order self-adjoint singularly perturbed boundary value problems. To demonstrate the competence of the method, we applied it on two model examples by taking different values for the perturbation parameter, ε, and mesh size, h. Numerical results obtained by the present method have been associated with numerical results obtained by the methods in [16, 17], and the results are summarized in Tables 1 and 4. Moreover, the maximum absolute errors decrease rapidly as the number of mesh points N increases. Further, as shown in Figs. 1 and 2, the proposed method approximates the exact solution very well for h ≥ ε, for which most of the current methods fail to give good results. To further verify the applicability of the planned method, graphs were plotted aimed at Examples 1 and 2 for exact solutions versus the numerical solutions obtained. As Figs. 1 and 2 indicate good agreement of the results, presenting exact as well as numerical solutions, which proves the reliability of the method. Also, Figs. 3 and 4 specify the effects of perturbation parameter and mesh sizes of the solution domain.
Table 1 Comparison of maximum absolute errors for Example 1
The behavior of exact and numerical solution for Example 1 at ε = 10−3 and N = 100
The physical behavior of the solution for Example 2 at N = 64 and ε = 10−3
Point-wise absolute errors of Example 1 at ε = 2−12 with different mesh size h
Point-wise absolute errors for Example 1 at N = 32 and different perturbation parameters
Further, the numerical results presented in this paper validate the improvement of the proposed method over some of the existing methods described in the literature. Both the theoretical and numerical error bounds have been established for the fourth and sixth-order methods. Hence, the Richardson extrapolation method accelerates fourth order into sixth order convergent as given in Table 2. The results in Tables 3 and 5 further confirmed that the computational rate of convergence and theoretical estimates are in agreement (Tables 4 and 5). Generally, the present method is consistent, is stable, and gives more accurate numerical solution for solving second-order self-adjoint singularly perturbed boundary value problems.
Table 2 Comparison between with and without Richardson methods of maximum absolute errors for Example 1
Table 3 Comparison of rate of convergence for Example 1
Table 5 Rate of convergence for Example 2
All data generated or analyzed during this study are included.
Delkhon, M, Delkhosh, M: Analytic solution of some self-adjoint equations by using variable change method and its application, J Appl Math, 2012, 1-7 (2012), http://dx.doi.org/ https://doi.org/10.1155/2012/180806.
Miller, HJJ, O'Riordan, E, Shishkin IG: Fitted numerical methods for singular perturbation problems, Error estimate in the maximum norm for linear problems in one and two dimensions, World Scientific. ISBN: 981-02-2462-1, (1996).
Suayip, Y., Niyazi, S.: Numerical solutions of singularly perturbed one-dimensional parabolic convection–diffusion problems by the Bessel collocation method. Appl Math Comput. 220, 305–315 (2013)
Roos, GH, Stynes, M, Tobiska, L: Robust numerical methods for singularly perturbed differential equations, Convection-diffusion-reaction and flow problems, Springer-Verlag Berlin Heidelberg, Second Edition, (2008).
Mishra, H.K., Kumar, M., Singh, P.: Initial-value Technique for self-adjoint singular perturbation boundary value problems. Comput Math Model. 20, 207–218 (2009)
Rashidinia, J., Mohammadi, R., Moatamedoshariati, S.H.: Quintic spline method for solution of singularly perturbed boundary value problems. Int J Comput Meth Eng Sci Mech. 11, 247–257 (2010)
Tirmizi, A.I., Haq, F.I., Islam, S.: Non-polynomial spline solution of singularly perturbed boundary value problems. Appl Math Comput. 196, 6–16 (2008)
Rashidinia, J., Mahmoodi, Z., Ghasemi, M.: Cubic spline solution of singularly perturbed boundary value problems with significant first derivatives. Appl Math Comput. 190, 1762–1766 (2005)
Kadalbajoo, MK, Kumar, D: Variable mesh finite difference method for self-adjoint singularly perturbed two-point boundary value problems, J Comput Math, (2010) doi: https://doi.org/10.4208/jcm.1003-m2809.
Kadalbajoo, M.K., Aggrwal, V.K.: Fitted mesh B-spline collocation method for solving self-adjoint singularly perturbed boundary value problems. Appl Math Comput. 161, 973–987 (2005)
Munyakazi, J.B.: Higher order numerical methods for singular perturbation problems, PhD thesis. University of the Western Cape, South Africa (2009)
Fasika, W.G., Gemechis, F.D., Tesfaye, A.B.: Tenth order compact finite difference method for solving singularly perturbed 1D reaction - diffusion equations. Int J Eng Appl Sci (IJEAS). 8, 15–24 (2016)
Fasika, W.G., Gemechis, F.D., Tesfaye, A.B.: Sixth-order compact finite difference method for singularly perturbed 1D reaction diffusion problems. J Taibah Univ Sci. 11, 302–308 (2017)
Fasika, W., Gemechis, F., Tesfaye, A.: Fourth order compact finite difference method for solving singularly perturbed 1D reaction diffusion equations with dirichlet boundary conditions. Momona Ethiop J Sci (MEJS). 8, 168–181 (2016)
Feyisa, E., Gemechis, F.: Higher order compact finite difference method for singularly perturbed one dimensional reaction diffusion problems. J Niger Math Soc. 36, 491–502 (2017)
Terefe, A., Gemechis, F., Tesfaye, A.: Fourth-order stable central difference method for self-adjoint singular perturbation problems. Ethiop J Sci Technol. 9, 53–68 (2016)
Yitbarek, Z., Gemechis, F., Tesfaye, A.: Sixth-order stable central difference method for self-adjoint singular perturbation problems. Ethiop J Educ Sc. 13, 23–41 (2017)
Zhilin, Li, Zhonghua Q, Tang, T: Numerical solution of differential equations, Introduction to finite difference and finite element methods, printed in the United Kingdom by Clays, (2018).
Kadalbajoo, M.K., Reddy, Y.N.: A non-asymptotic method for general singular perturbation problems. J Optimization Appl. 55, 256–269 (1986)
Smith, GD: Numerical solution of partial differential equations, Finite difference methods, Third edition, Oxford University Pres, New York (1985).
The authors wish to express their thanks to the authors of literatures for the provided scientific aspects and idea for this work. Also, we request to express great thanks to reviewers for the constructive comments and inputs given to improve the quality of our work.
Department of Mathematics, Jimma University, P. O. Box 378, Jimma, Ethiopia
Muslima Kedir Siraj, Gemechis File Duressa & Tesfaye Aga Bullo
Muslima Kedir Siraj
Gemechis File Duressa
Tesfaye Aga Bullo
MKS raised the first idea on the stated title and wrote, analyzed and interpreted the data. GFD modified the title so as to make it workable, edited, analyzed and interpreted the data again in the well-organized form. TAB accomplished the numerical computation of the solution for the governing differential equations of the present problem using MATLAB software and display results in tables and graphs. All authors read and approved the final manuscript.
Correspondence to Tesfaye Aga Bullo.
Fourth-order stable central difference with Richardson extrapolation method has been formulated for solving second-order self-adjoint singularly perturbed boundary value problems using the study design of both documentary review and numerical experimental using MATLAB R2013a software which gives more accurate numerical solution with the corresponding sixth-order convergent.
Siraj, M.K., Duressa, G.F. & Bullo, T.A. Fourth-order stable central difference with Richardson extrapolation method for second-order self-adjoint singularly perturbed boundary value problems. J Egypt Math Soc 27, 50 (2019). https://doi.org/10.1186/s42787-019-0047-4
Singular perturbation
self-adjoint problem
boundary value problem
Mathematics subject classification
65L10, 65L11, 65L12, 65B05, 65Y04