A NSFD method for the singularly perturbed Burgers-Huxley equation

This article focuses on a numerical solution of the singularly perturbed Burgers-Huxley equation. The simultaneous presence of a singular perturbation parameter and the nonlinearity raise the challenge of finding a reliable and efficient numerical solution for this equation via the classical numerical methods. To overcome this challenge, a nonstandard finite difference (NSFD) scheme is developed in the following manner. The time variable is discretized using the backward Euler method. This gives rise to a system of nonlinear ordinary differential equations which are then dealt with using the concept of nonlocal approximation. Through a rigorous error analysis, the proposed scheme has been shown to be parameter-uniform convergent. Simulations conducted on two numerical examples confirm the theoretical result. A comparison with other methods in terms of accuracy and computational cost reveals the superiority of the proposed scheme.

The field of numerical methods for singularly perturbed problems has flourished significantly over the last couple of decades. The studies focused on the design and analysis of numerical methods that were parameter-uniformly convergent. In other words, research was mostly motivated by the challenge posed by the presence of the singular perturbation parameter. A vast majority of researchers discussed linear problems. Nonlinear singularly perturbed equations, and in particular Equation (1.1), received little attention from scholars in this research area. Nonlinearity constitutes an another layer of difficulties of singularly perturbed problems.
In general, the solution of the singularly perturbed Burgers-Huxley problem (1.1) with arbitrary initial and boundary conditions (IBCs) cannot be expressed in terms of a finite number of elementary functions [2,3]. Thus, several scholars seek approximate analytic solutions using .
/fams. . different transformation techniques [4][5][6]. However, these techniques still hold for only some specific parameters and initial boundary conditions. Hence, for general cases, there is a need to develop numerical methods to obtain an approximate solution of the problem. For the last two decades, various numerical methods have been studied to solve the Burgers-Huxley equations of type (1.1) for ε = 1 and specific IBCs in the framework of NSFD methods. For example, A.R. Appadu et al. [7] developed two novel nonstandard finite difference schemes, and explicit exponential finite difference method and a fully implicit exponential finite difference method.
For the singularly perturbed case, that is when 0 < ε ≪ 1, various parameter-uniformly convergent methods have been developed to solve (1.1) [1,2,8,9]. It is to be noted that all these works used the Newton's quasilinearization to deal with the nonlinearity.
In this article, a nonstandard finite difference method for the singularly perturbed case is proposed. To deal with the nonlinearity, rather than using quasilinearization, for the first time, nonlinear terms are approximated in a nonlocal manner following one of Micken's rules [10]. The resulting method preserves the properties of the continuous solution and provides reliable numerical results. The method is proved to be first order parameter-uniform convergent in time and space.
The rest of this article is structured in the following manner. In the next section, we study a priori estimates of the analytic solution to the problem. Section 3 is about the proposed NSFD method and its parameter-uniform convergence. Section 4 deals with the implementation of the method to confirm the theoretical results and compare with other methods. Section 5 concludes the present work and provides a direction for future work.
. Some analytical results: A priori estimates Throughout this article, we assume that u 0 , h 1 , and h 2 are sufficiently smooth functions and that a and b satisfy Under these assumptions the problem (1.1) has a unique solution u(x, t) ∈ C 2,1 (D) which exhibits a boundary layer near x = 1 [1].
In the rest of the paper, C denotes a generic constant independent of the parameters and mesh sizes, and . D denotes the maximum norm over D. For the solution u(x, t) of (1.1), we have the following bounds.

Lemma 2.1. [Uniform stability estimate for continuous problem]
Let u(x, t) be the exact solution of (1.1) onD, then Proof: The proof can be seen in [9]. Remark 2.2. From our assumption in (2.1) and Lemma 2.1, one can observe that the solution u(x, t) satisfies . Proposed scheme and its convergence analysis This section is dedicated to the construction of a new scheme and to the analysis of its parameter-uniform convergence. The first step is to discretize the solution domainD and provide the definition of the NSFD methods given in [11] (a revised version of [10]).
Let  1. The classical denominator h or h 2 of the discrete first or second order derivative is replaced by a nonnegative function ψ such that For example, denominator functions that satisfy the above conditions are and so on. and so on.

Definition 3.2.
Assume that the solution of (1.1) satisfies some property P. The difference equation of (1.1) in u n m is called (qualitatively) stable with respect to P if, for any values of the mesh sizes t and h, solution of the difference equation replicates the property P. The remaining rules are expressed in terms of definition 3.2. For example, the schemes under consideration in this paper is qualitatively stable with respect to the order of the differential equation, they do satisfy positivity and uniform boundedness.
To construct the scheme, first we semidiscretize the problem (1.1) in time direction and then in space direction.

. . Semidiscrete scheme
Denote u n (x) as the approximation of u(x, t n ) at time level t n , 0 ≤ n ≤ N. Now, we apply backward Euler finite difference method to discretize the continuous problem (1.1) in the temporal direction and obtain the following semidiscrete scheme: Here the nonlinearity a(u)u x and b(u)u approximated in a nonlocal way as and thus this approximation satisfies the condition (2) in definition 3.1.

. . Convergence analysis for the semidiscrete scheme
The local truncation error e n+1 of the semidiscrete scheme (3.3) is given by e n+1 = u(x, t n+1 ) −û n+1 (x), whereû n+1 (x) is the solution of (3.4) That is,û n+1 (x) is the solution obtained after one step of semidiscrete scheme (3.3) by taking the exact value u(x, t n ) instead of u n (x) as the starting data.
In order to analyze the uniform convergence of the solution u n (x) of (3.3) to the exact solution u(x, t n ), we will perform the stability analysis and establish the consistency result. First, let us consider the semidiscrete maximum principle for the operator I + tL N ε .
From the given hypothesis and second derivative Assumption (2.1) leads to (I + tL N ε ) n+1 (x * ) < 0, which contradicts the given assumption in Lemma 3.4 and thus n+1 (x) ≥ 0 for all x ∈¯ . This maximum principle leads to the following stability result

Lemma 3.5. (Local error estimate)
Estimate of the local error e n+1 is given by The global error E n+1 associated to the semidiscrete scheme (3.3) at (n + 1)-th time level is given by Using the local error estimate (Lemma 3.5) and triangular inequality the following global error estimate follows.

Theorem 3.6. (Global error estimate)
The global error E n+1 of the time discretisation at the n + 1 time step satisfies Therefore, the semidiscrete scheme (3.3) is a first order uniformly convergent in the time direction.
The following lemma provides the asymptotic estimates of the exact solution u n+1 of (3.3) and its derivatives. These estimates will be used in the convergence analysis of the fully discrete scheme (3.8) Proof: See the proof of Lemma 4.1 in [12] . . The spatial discretisation In this subsection, we totally discretise the semidiscrete scheme (3.3) on a uniform mesh¯ M . Let the approximation of u n (x) at x m be denoted by u n m , 0 ≤ m ≤ M. Similarly, let the approximations of a(u n m ) and b(u n m ) be denoted respectively by a n m and b n m . Using the upwind finite difference scheme and nonstandard finite difference methodology of Mickens [13], the semidiscrete problem (3.3) can be fully-discretised as      Equation (3.10) leads to the tridiagonal system which can be solved with Thomas Algorithm [15,16].

. . Convergence analysis of the space discretisation
The difference operator I + tL M,N ε of (3.9) satisfies the following maximum principle. Hence (3.9) has a unique discrete solution u n m .  This leads immediately to the following stability result, analogous to the continuous result.  Proof: Define Proof: The truncation error of the complete discretization (3.9) is given by (3.15) where Applying absolute values and using triangle inequalities in (3.15) leads to From the use the Taylor series expansions of some terms in the above equation, on obtains h 2 (ψ n m ) 2 − 1 = − a n (η 1 )h 2ε + ha(η 1 ) , (3.17) for some η i such that x m−1 ≤ η i ≤ x m+1 , i = 1, 2, 3, 4.  (3.20) in (3.16) and from the boundedness of a n (x) in (3.8), one has m ) ≤ ε t a n (η 1 )h 2ε + ha n (η 1 ) +ε t h 2 16 (3.21) and this gives Applying Lemma 7 of [17] gives (3.24) The temporal and spatial error estimates in Theorem 3.6 and 3.12, respectively, give the the following main result of this paper.
Therefore, the presented discrete scheme is ε-uniform convergent of order one both in time and space.

. Numerical implementation
In this section, two test examples are provided to demonstrate the efficiency and applicability of the proposed numerical method. The first example is taken from a recent article [7] for ε = 1, and we modified it by multiplying the highest derivative term by ε to make the problem singularly perturbed. It is the first time to consider this example for the singularly perturbed case. The second example is taken from [2] for different initial and boundary conditions to satisfy our assumption (2.1). All the computations are carried out by Intel Coreå i5-4210M CPU @2.60GHz × 4 with MATLAB 2017.
For ε = 10 −4 , α = β = 1 and γ = 0.5, we plotted the numerical solution of Example 4.1 using the proposed method in Figure 1. As we see, the method resolves the boundary layer at x = 1. Thus, our proposed method replicates the property of the continuous solution.
Since the exact solution of this problem is not known, for each ε, the maximum pointwise error is calculated using double mesh principle [18] given by where u n;N m;M and u 2n;2N 2m;2M are approximate solutions to problem (3.9) on D M,N and D 2M,2N , respectively. The corresponding order of convergence is given by In Tables 1, 2, we compute the maximum absolute errors and the corresponding order of convergences using the proposed method for α = β = 1, γ = 0.5 and various values of ε with fixed values of M and N, respectively. From these results, we can observe that the method is ε-uniform convergent of order one both in time and space. This confirms the theoretical error results. Table 3 compares absolute errors at some values of x and t using the proposed method and methods in [7,19] and the nonlinear nonstandard finite difference (NNSFD) method given by where the denominator functions The nonlinear system (4.1) can be easily solved using Newton's method (see in [7]). The solution at the previous time-step is taken as the initial estimate. The methods considered in [7] and [19] are NSFD methods and variational iteration method (VIM), respectively. From the table, we observe that our proposed method gives more accurate results.

. Conclusion
In this paper, a nonstandard finite difference method for a singularly perturbed Burgers-Huxley equation has been developed. First, the backward-Euler scheme was applied to discretize the problem (1.1) with respect to time derivative and the upwind nonstandard finite difference scheme on uniform mesh to approximate the spatial derivative. Then, the presented method was proved to be first-order convergent in both the spatial and temporal variables. Numerical results are given in Figures 1, 2 Tables 1-6 for two test examples to confirm the theoretical results and to compare with recent results. It has been observed from these figures and tables that the numerical results are in agreement with the theoretical findings. For ε = 1, in Table 3, comparisons with the NNSFD and the VIM of [19] and the NSFD method of [7] reveal that the proposed NSFD method gives more accurate results. In addition, the present method is also applicable when 0 < ε < 1. Thus, in all cases, the present method produces more accurate results than the existing schemes.

and in
For future work, one may proceed with a higher order scheme for the problem under consideration for example by using the Crank-Nicolson method or those presented in [20,21]. We are currently working in this direction. Also to note is that higher order parameteruniform numerical methods for Burgers-Huxley equations are quasi-absent. The design of direct higher order parameteruniform convergent methods is thus an interesting direction to explore.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions
ED constructed the scheme and implemented the MATLAB code under JM's supervision. All authors participated in the analysis of the theoretical results and in the writing of the manuscript.