Uniformly convergent computational method for singularly perturbed unsteady burger-huxley equation

Graphical abstract

Many researchers have made efforts to construct analytical and numerical methods for Burger equations, for instance, [6,[8][9][10][11]16,17,[19][20][21]25,28,29,32] and references therein. A Burger-Huxley equation in which the highest order derivative is multiplied by a small perturbation parameter ε (0 < ε << 1) is termed as the singularly perturbed Burger-Huxley equation. Due to the presence of perturbation parameter ε the solution reveals boundary/ sharp interior layer(s), and it is not easy to find a stable numerical approximation.
The methods suggested by [6,[8][9][10][11]16,17,[19][20][21]25,28,29,32] and other classical numerical methods on a uniform mesh fail to approximate the singularly perturbed Burger-Huxley equation. They require an unacceptably large number of mesh points to sustain the approximation because the mesh width depends on the ε. This limitation of the conventional numerical methods has encouraged researchers to develop robust numerical techniques that perform well enough independent of the ε. [14] presented a uniformly convergent finite difference method (FDM) for the problem (1.1) .
[7] developed a numerical method that comprises of an implicit-Euler method to discretize in temporal direction on a uniform mesh and a monotone hybrid finite difference operator to discretize the spatial variable with a piecewise uniform Shishkin mesh.
A robust numerical method that consists of backward-Euler scheme on a uniform mesh to approximate time derivative and upwind FDM on an adaptive nonuniform grid is used for space derivative of Eq. (1.1) is suggested by [23] . [12,13] considered a similar problem as in [23] and proposed a parameter uniform numerical method based on fitted operator techniques. Very recently, [3] developed fitted operator finite difference method for singularly perturbed Burgers' initialboundary value problem.
Nevertheless, the numerical treatment of problem (1.1) has seen little development. The purpose of this work is to formulate and analyze an ε−uniformly convergent numerical method for solving the problem (1.1) . The novelty of the presented method, unlike Shishkin mesh developed by [30] and Bakhvalov mesh developed by [2] , does not require a priori information about the location and width of the boundary layer. The methods developed by [7,[12][13][14]23] for Eq. (1.1) are based on finite difference method. In comparison with the finite difference methods, spline solution has its own advantages. For instance, once the solution has been computed, the information required for spline interpolation between mesh points is available [15] . Also, splines are a simpler and more practical way to solve boundary-value problems than finite difference methods [18] . This is particularly important when the solution of the boundary-value problem is required at various locations in the interval [ a, b] . This provides the motivation for our work on using fitted exponential spline method for solving the problem under consideration.
A priori estimates for the solution of the continuous problem Proof. See [7] .

Quasi-linearization
To linearize the semi-linear term of Eq. (1.1) , we choose the reasonable initial approximation for the function (z, t ) in the term ζ (z, t, (z, t ) , z (z, t )) , and denote it as 0 (z, t ) that satisfy both initial and boundary conditions which is obtained by separation of variables method of the homogeneous part of the problem under consideration and is given by [13] : Thus, the nonlinear term ζ (z, t, (z, t ) , z (z, t )) of Eq. (3.1) can be linearized by applying the Newton-Raphson-Kantorovich approximation approach as

Temporal semi-discretization
To discretize the time variable for Eq. (3.4) , we use the implicit Euler method with uniform time (3.5) exhibits boundary located at z = 1 and admits a unique solution. Clearly, the operator 1 + τ £ M ε satisfies the maximum principle, which ensures the stability of the semi-discrete (3.5) .
where C is a positive constant independent of ε and t.
Proof. The use of Taylor's series expansion for ˆ k +1 (z, t j+1 ) provides where ˆ k +1 (z, t j+1 ) is the solution of the boundary value problem (3.5) .
Hence, applying the maximum principle on the operator gives

Lemma 3.2 (Global Error Estimate) . Under the hypothesis of Lemma 3.1 , the global error estimate in the temporal direction is given by
Proof. The use of Lemma 3.1 yields where c 1 , C are positive constants independent of ε and τ .

Lemma 3.3. The bound on the derivative of the solution ˆ k +1 (z) of the Eq. (3.5) is given by
Proof. See [7] .

Spatial Semi-discretization
The space domain [0,1] is divided into N equidistant mesh with the spatial step size such that Let us rewrite Eq. (3.5) for i = 0 , 1 , · · · , N, and j = 0 , 1 · · · , M − 1 as To solve Eq. (3.9) , we use an exponential cubic spline method as discussed below. Let ˆ j+1 i be an where, A i , B i , G i , and D i are constants and ψ = 0 is a free parameter that will be used to enhance the accuracy of the method. To find the unknown coefficients in Eq. (3.10) , let us denote: (3.11) Differentiating twice both sides of Eq. (3.10) with respect to z yields (3.12) Using the relation in Eq. (3.11) into (3.12) at z = z i , and it yields Again substituting Eq. (3.11) into (3.12) at z = z i +1 , and it gives (3.14) From Eqs. (3.13) and (3.14) , Using the continuity condition of the first derivative at where Equation (3.9) can be rewritten as where we approximate ˆ j+1 Now, introducing a fitting factor σ ( ρ) , which is used to handle the effect of the ε on solution behavior, into the above Eq. (3.21) provides where ρ = ε .
Taking the limit of Eq. (3.22) as → 0 : (3.23) For problems with a layer at the right end of the interval, using the theory of singular perturbations, the solution of Eq. (3.9) takes the form [27] where j+1 0 (z) is the solution of the reduced problem Considering Taylor's series expansion for j+1 (z) about the point z = 1 and restricting to their first terms, Eq. (3.24) becomes Plugging the above equations into Eq. (3.23) , we obtain the required fitting factor (3.26) Finally, using Eqs. (3.22) and Eq. (3.26) : For sufficiently small mesh sizes the above matrix is non-singular and the matrix are diagonally dominant). Hence, by [26] , the matrix ϑ is M-matrix and have an inverse. Therefore, the system of equations can be solved by matrix inverse with the given boundary conditions. This lemma gives guarantee for the existence of unique discrete solution.

Convergence analysis
On the boundary points: Now, on the discretized spatial domain ϒ N : By Lemma 4.1 , we obtain Π j+1 i ± ≥ 0 , 0 ≤ i ≤ 1 . Hence, the required bound is obtained.

Lemma 4.3. The local truncation error in space semi-discretization of the discrete problem (3.27) is given as
where C is a constant independent of ε and .

(4.2)
Considering the corresponding exact solution to Eq. (4.2) : (4.5) Hence, for the choice of ω 1 + ω 2 = 1 / 2 , we obtain The matrix representation of Eq. (4.6) is Following [31] , it can be shown that where C is a constant, independent of and ε.  Let (z, t ) be the solution of problem (3.4) at each grid point z i , t j+1 and ˆ j+1 i be its approximate solution obtained by the proposed scheme given in Eq. (3.27) . Then, the error estimate for the fully discrete method is given by

Numerical examples, results and discussion
Several model examples have been presented to illustrate the efficiency of the proposed method. As the exact solutions of the considered examples are not known, we calculate the maximum absolute error for each ε given in [4] by: where z i , t j+1 denotes the numerical solution obtained in ϒ N,M with N mesh intervals in the spatial direction and M mesh intervals in the temporal direction. We determine the corresponding order of convergence for each ε by: Now, for each N and τ , the ε−uniform maximum point wise error is calculated using: and the corresponding ε−uniform order of convergence is calculated by: The computed E N,τ ε , E N,τ and R N,τ of the considered problems using the proposed scheme are depicted in Tables 1-7 . From these tables, one can observes that the developed scheme converges independent of the perturbation parameter with order of convergence one. Concisely, the numerical results show that the proposed method provides better results than results in [14,23] . In Fig. 1 , the 3D view of the solution of Examples 5.1 and 5.2 are given for visualizing the boundary layer at N = 64 , M = 40 and ε = 2 −18 . In Fig. 2 , the effect of ε on the solution profile is shown. As one observes, as ε → 0 strong boundary layer is created near z = 1 . In Fig. 3 , the effect of the time step on the solution profile is shown. This figure shows the existence of the boundary layer at z = 1 with time variable t → 1 . Fig. 4 shows the uniform convergence of the proposed scheme in log-log scale plot.     [14] 4.1129e-2 2.2723e-2 1.2008e-2 6.179e-3 3.136e-3

Conclusion
In this paper, uniformly convergent computational scheme is developed, analyzed, and implemented for solving singularly perturbed non-linear Burgers' problem which has a boundary layer. The developed scheme constitutes the implicit Euler in the time discretization and the fitted exponential cubic spline method in space discretization. It is shown that the developed scheme converges ε−uniformly which is first-order accurate in time and second-order accurate in space.
Four test examples are considered to validate the theoretical findings. The results in the tables and figures indicate the parameter uniform convergence of the scheme with the rate of convergence of one. Also, the experiment shows that the proposed scheme provides better numerical accuracy compared to schemes in [3,14,23] . From the results presented, one can conclude that the proposed scheme is an efficient robust numerical scheme to approximate the solution of the singularly perturbed unsteady non-linear Burgers' equation. The scheme can also be extendable for solving higher order families of singularly perturbed partial differential problems and singularly perturbed problems with degenerate coefficients. The obtained numerical results and graphs are plotted by using the MATLABR2013a software package.

Funding statement
This research did not receive any specific grant from funding agencies.

Declaration of Competing Interest
The authors declare that they have no conflicts of interest.