## ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 04 January 2023
Sec. Numerical Analysis and Scientific Computation
Volume 8 - 2022 | https://doi.org/10.3389/fams.2022.1061245

# Numerical treatment of singularly perturbed unsteady Burger-Huxley equation

Imiru Takele Daba1* Gemechis File Duressa2
• 1Department of Mathematics, Dilla University, Dilla, Ethiopia
• 2Department of Mathematics, Jimma University, Jimma, Ethiopia

This article deals with the numerical treatment of a singularly perturbed unsteady Burger-Huxley equation. The equation is linearized using the Newton-Raphson-Kantorovich approximation method. The resulting linear equation is discretized using the implicit Euler method and an exponential spline method for time and space variables, respectively. Richardson's extrapolation technique is employed to increase the accuracy in the temporal direction. The stability and uniform convergence of the proposed scheme are investigated. The scheme is shown uniformly convergent with the order of convergence O(τ + ℓ2) and O2 + ℓ2) before and after Richardson extrapolation, respectively. Several test examples are considered to validate the applicability and efficiency of the scheme. It is observed that the proposed scheme provides more accurate results than the methods available in the literature.

## 1. Introduction

Singular perturbed non-linear parabolic partial differential equations have many applications in real life scenarios. A good example of singularly perturbed nonlinear differential equations is the mathematical model for the prototype model used to describe the interaction between non-linear convection effects, reaction mechanisms, and diffusion transport. This equation has many intriguing phenomena such as bursting oscillation [1], interspike [2], population genetics [3], bifurcation, and chaos [4]. Several membrane models based on the dynamics of potassium and sodium ion fluxes were found in Lewis et al. [5].

Many researchers have made efforts to construct analytical and numerical methods for Burger equations, for instance, Ismail et al. [6], Javidi et al. [7], Estevez [8], Krisnangkura et al. [9], Satsuma et al. [10], Wang et al. [11], Hashim et al. [12, 13], Khattak [14], Krisnangkur et al. [9], Mohammadi [15], Sari et al. [16], Appadu et al. [17, 18], and Appadu and Tijani [19], and references therein. A Burger-Huxley equation in which the highest order derivative is multiplied by a perturbation parameter ε (0 < ε << 1) is termed the singularly perturbed Burger-Huxley equation (SPBHE). Due to the presence of the perturbation parameter and the nonlinearity in the problem, finding their solutions are tedious task. For instance, due to the presence of the perturbation parameter, the solution exposes boundary/shrill interior layer(s), and it is not easy to find a stable numerical approximation. The methods suggested in Ismail et al. [6], Javidi and Golbabai [7], Estevez [8], Krisnangkura et al. [9], Satsuma et al. [10], Wang et al. [11], Hashim et al. [12, 13], Khattak [14], Mohammadi [15], and Sari and Gürarslan [16] and other standard numerical methods on a uniform mesh is insufficient to approximate the SPBHE. Thus, to resolve this drawback, the fitting approaches are competitive computational techniques to overcome the limitations of classical numerical methods. Kaushik and Sharma [20] presented a uniformly convergent finite difference method (FDM) for the problem Equation (1). Gupta and Kadalbajoo [21] developed a numerical method that comprises an implicit-Euler method and a monotone hybrid finite difference operator with a Shishkin mesh type for the problem Equation (1). A robust numerical method that consists of a backward-Euler method and an upwind FDM on an adaptive nonuniform grid for Equation (1) is suggested by Liu et al. [22]. Kabato and Duressa [23] and Jima et al. [24] considered a similar problem as in Liu et al. [22] and proposed a parameter uniform numerical method based on fitted operator techniques. Kehinde et al. [25] developed a fitted operator FDM for the two-dimensional semilinear singularly perturbed convection-diffusion problem.

However, there has been little development in the numerical treatment of the problem Equation (1). This work aims to develop and analyze an ε, which is a uniformly convergent numerical method for solving the problem Equation (1). First, the nonlinear terms are linearized by using the Newton-Raphson-Kantorovich approximation method. Then, we approximate the resulting singularly perturbed parabolic convection-diffusion equations using the implicit Euler method for the time derivative and fitted the exponential spline method for the space derivative on a uniform step length. Finally, to enhance the order of convergence in the time variable, the Richardson extrapolation technique is applied to the temporal variable. The novelty of the presented method, unlike Shishkin and Bakhvalov mesh types, does not require a priori information about the location and width of the boundary layer.

## 2. Continuous problem

On the domain ϒ = χz × χt = (0, 1) × (0, T], we consider singularly perturbed Burger-Huxley equation of the form

Where 0 < ε ≪ 1 is a singular perturbation parameter and α ≥ 1, θ ≥ 0, and λ ∈ (0, 1). The functions Θ0(t), Θ1(t), and $\Im$0(z) are assumed to be sufficiently smooth, bounded, and independent of ε.

Lemma 2.1 (Maximum principle). If $\Im \in {C}^{2,1}\left(\overline{\Upsilon }\right)$, $\Im$|∂ϒ ≥ 0, and £εℑ|ϒ ≥ 0, then $\Im {|}_{\overline{\Upsilon }}\ge 0$.

Proof. See Gupta and Kadalbajoo [21].

An immediate consequence of Lemma 2.1 for the solution of problem Equation (1) resulted in Lemma 2.2.

Lemma 2.2 (Stability estimate). Let $\Im$(z, t) be the solution of problem Equation (1), then we have,

$||ℑ||ϒ¯≤T||ℑ0||Γi+||ℑ||∂ϒ.$

Proof. See Gupta and Kadalbajoo [21].

## 3. Formulation of the numerical scheme

### 3.1. Quasi-linearization

The Equation (1) is rewritten as

Where $\zeta \left(z,t,\Im \left(z,t\right),{\Im }_{z}\left(z,t\right)\right)=-\alpha \Im \frac{\partial \Im }{\partial z}+\theta \left(1-\Im \right)\left(\Im -\lambda \right)$ is the non-linear function of z, t, $\Im$(z, t), $\Im$z(z, t).

To linearize the semi-linear term of Equation (1), we choose the reasonable initial approximation for the function $\Im$(z, t) in the term ζ(z, t, $\Im$(z, t), $\Im$z(z, t)) and denote it as $\Im$0(z, t) that satisfy both initial and boundary conditions which are obtained by the separation of variables method of the homogeneous part of the problem under consideration and is given by Kabeto and Duressa [23].

Thus, the nonlinear term ζ(z, t, $\Im$(z, t), $\Im$z(z, t)) of Equation (2) can be linearized by applying the Newton-Raphson- Kantorovich approximation approach as

Where ${\left\{{\Im }^{\left(k\right)}\right\}}_{k=0}^{\infty }$ are a sequence of approximate solutions of ζ(z, t, $\Im$(z, t), $\Im$z(z, t)).

Substituting Equation (4) into Equation (2), we have

where

### 3.2. Temporal semi-discretization

We define the uniform mesh for the domain χt as ${\chi }_{t}={\Upsilon }_{\tau }^{M}=\left\{{t}_{j}=j\tau ,\forall j=0,2,3,\cdots \phantom{\rule{0.3em}{0ex}},M,\tau =1/M\right\}$ and we approximated the temporal derivative term of Equation (5) using implicit Euler method, which gives

Since $\varpi \left(z,{t}_{j+1}\right)\ge {\varpi }^{*}>0$ and $q\left(z,{t}_{j+1}\right)\ge {q}^{*}>0,z\in {\overline{\chi }}_{z}$ Equation (6) exhibits boundary located at z = 1 and admits a unique solution.

Lemma 3.1 (Local Error Estimate). In the solution of Equation (6), the local truncation error estimate is bounded by

$||LTEj+1||∞≤C*(τ)2,$

Where C* is a positive constant independent of ε and τ.

Proof. On applying Taylor's series expansion on ${\stackrel{^}{\Im }}^{k+1}\left(z,{t}_{j+1}\right)$ yields

$ℑ^k+1(z,tj+1)=ℑ^k+1(z,tj)+τℑ^tk+1(z,tj)+O((τ)2).$

This implies

Substituting Equation (5) into Equation (7) we have

Subtracting Equation (6) from Equation (8), and the local truncation error $LT{E}_{j+1}={\stackrel{^}{\Im }}^{k+1}\left(z,{t}_{j+1}\right)-{\stackrel{^}{\Im }}^{k+1}\left(z,{t}_{j+1}\right)$ at (j + 1)th is the solution of a BVP

Hence, applying the maximum principle on the operator gives

$||LTEj+1||∞≤C*(τ)2.$

Lemma 3.2. [Global Error Estimate (GEE)] The GEE in the temporal direction of Equation (6) is given by

$||Ej||∞≤C(τ),∀j≤T/τ.$

Proof. Using Lemma 3.1, we obtain the following GEE at (j)th time step

Where c1, C are positive constants independent of ε and τ.

Lemma 3.3. The solution ${\stackrel{^}{\Im }}^{k+1}\left(z\right)$ of the Equation (6) is bounded by

Proof. See Gupta et al. [21].

#### 3.2.1. Spatial semi-discretization

We define the uniform mesh for the space domain [0, 1] as

$ϒℓN={zi=iℓ,∀i=1,2,3,⋯,N,z0=0,zN=1,ℓ=1/N}.$

Equation (6) can be rewritten as

Where $r\left(z,{t}_{j+1}\right)=\left(q\left(z,{t}_{j+1}\right)+\frac{1}{\tau }\right)$ and $\phi \left(z,{t}_{j+1}\right)=v\left(z,{t}_{j+1}\right)+\frac{{\stackrel{^}{\Im }}^{k+1}\left(z,{t}_{j}\right)}{\tau }$.

For spatial discretization of Equation (10), we use an exponential cubic spline method.

The approximation solution ${\stackrel{^}{\Im }}_{i}^{j+1}$ of Equation (10) obtained by the segment S(z) passing through the points $\left({z}_{i},{\stackrel{^}{\Im }}_{i}^{j+1}\right)$ and $\left({z}_{i+1},{\stackrel{^}{\Im }}_{i+1}^{j+1}\right)$. Each mixed spline segment ${S}_{△}^{j+1}\left(z\right)$ has the form [26].

Where Ai, Bi, Ci, and Di are constants and ψ ≠ 0 is a free parameter that is used to enhance the accuracy of the method. To find the values of Ai, Bi, Ci, and Di in Equation (11), we denote

Differentiating twice both sides of Equation (11), we get

Using the relation in Equation (12) into Equation (13) at z = zi, we get

$d2S△j+1(zi)dz2=Ξi=ψ2(Ai+Bi)$

and it yields

Where ϱ = ψℓ and i = 0, 1, 2, ⋯ , N.

Again substituting Equation (12) into Equation (13) at z = zi+1, we obtain

$d2S△j+1(zi+1)dz2=Ξi+1=ψ2(Aieϱ+Bie-ϱ)$

and it gives

From Equations (14) and (15), we have

Substituting Equation (12) into Equation (11) at z = zi and z = zi+1, we obtain

which implies

The first derivative continuity condition at z = zi, that is S△−1(zi) = S(zi) gives

where

Now from Equation (10), we have

Where we approximate ${\stackrel{^}{\Im }}_{z}^{j+1}\left({z}_{i-1}\right),{\stackrel{^}{\Im }}_{z}^{j+1}\left({z}_{i}\right)$ and ${\stackrel{^}{\Im }}_{z}^{j+1}\left({z}_{i+1}\right)$ as

Substituting Equation (20) into (19), we have

Substituting Equation (21) into Equation (18) and rearranging, we get

To grip the effect of the perturbation parameter ε, we multiply the perturbation parameter of Equation (22) by fitting factor σ(ρ), we obtain

Where $\rho =\frac{\ell }{\epsilon }$.

Evaluating the limit of Equation (23) as ℓ → 0

When the boundary layer is on the right side of the domain, from the theory of singular perturbation [27], the solution of Equation (10) is of the form

Where ${\Im }_{0}^{j+1}\left(z\right)$ is the solution of the reduced problem

Using Taylor's series expansion for ϖj+1(z) about the point z = 1 and restricting to their first terms, Equation (25) becomes

Equation (26) at zi = iℓ and as ℓ → 0 becomes

Plugging the above equations into Equation (24), we get

Finally, from Equations (23) and (27), we get

where

${£εN,Mℑ^ij+1=ϑi-ℑ^i-1j+1+ϑicℑ^ij+1+ϑi+ℑ^i+1j+1,ϑi-=-σ(ρ)ρ-3ω1ϖi-1j+12-ω2ϖij+1+ω1ϖi+1j+12+ℓω1ri-1j+1,ϑic=2σ(ρ)ρ+2ω1ϖi-1j+1-2ω1ϖi+1j+1+2ℓω2rij+1,ϑi+=--σ(ρ)ρ-ω1ϖi-1j+12+ω2ϖij+1+3ω1ϖi-1j+12+ℓω1ri+1j+1,Hij+1=ℓ(ω1φi-1j+1+ω2φij+1+ω1φi+1j+1).$

For sufficiently small mesh sizes, the above matrix is non-singular and $|{\vartheta }_{i}^{c}|\ge |{\vartheta }_{i}^{-}|+|{\vartheta }_{i}^{+}|$ (i.e., the matrix are diagonally dominant). Hence, by Nichols [28], the matrix ϑ is M-matrix and has an inverse. Therefore, the system of equations can be solved by matrix inverse with the given boundary conditions.

## 4. Convergence analysis

Lemma 4.1 (Discrete Maximum Principle). If the discrete function ${\Im }_{i}^{j+1}$ satisfies ${\Im }_{i}^{j+1}\ge 0$, on ∂ϒ. Then ${£}_{\epsilon }^{N,M}{\Im }_{i}^{j+1}\ge 0$ on ϒN,M implies that ${\Im }_{i}^{j+1}\ge 0$ at each point of ${\overline{\Upsilon }}^{N,M}.$

This lemma gives assurance for the presence of a unique discrete solution.

Lemma 4.2 (Discrete Uniform Stability). The solution ${\stackrel{^}{\Im }}_{i}^{j+1}$ of the discrete problem (28) at (j + 1)th time level and $\eta ={min}_{0\le i\le N}\left\{{r}_{i}^{j+1}\right\}$, where η is some positive constant that satisfies

$||ℑ^ij+1||≤||£εN,Mℑ^ij+1||η+max{|ℑ^0j+1|,|ℑ^Nj+1|}.$

Proof. We define barrier functions ${\left({\Pi }_{i}^{j+1}\right)}^{±}$ as

$(Πij+1)±=Z±ℑ^ij+1,$

Where $Z=\frac{||{£}_{\epsilon }^{N,M}{\stackrel{^}{\Im }}_{i}^{j+1}||}{\eta }+max\left\{|{\stackrel{^}{\Im }}_{0}^{j+1}|,|{\stackrel{^}{\Im }}_{N}^{j+1}|\right\}$.

On the boundary points, we obtain

Now, on the discretized spatial domain ${\Upsilon }_{\ell }^{N}$, we have

$£εN,M(Πij+1)±=£εN,M(Z±ℑ^ij+1)=(-σ(ρ)ρ-3ω1ϖi-1j+12-ω2ϖij+1+ω1ϖi+1j+12+ℓω1ri-1j+1)(Z±ℑ^i-1j+1)+(2σ(ρ)ρ+2ω1ϖi-1j+1-2ω1ϖi+1j+1+2ℓω2rij+1)(Z±ℑ^ij+1)+(--σ(ρ)ρ-ω1ϖi-1j+12+ω2ϖij+1+3ω1ϖi-1j+12+ℓω1ri+1j+1)(Z±ℑ^i+1j+1),$

By Lemma 4.1, we obtain ${\left({\Pi }_{i}^{j+1}\right)}^{±}\ge 0,0\le i\le 1$. Hence, the required bound is obtained.

Lemma 4.3. The local truncation error in space semi-discretization of the discrete problem (28) is given as

$||ℑ^j+1(zi)-ℑ^ij+1||≤Cℓ2,$

Where C is a positive constant independent of ε and ℓ.

Proof. From the truncation error of Equation (20), we have

Where zi−1 < ξ < zi+1.

Substituting

into Equation (18), we get

Considering the corresponding exact solution to Equation (30), we have

Subtracting Equation (30) from Equation (31) and denoting ${e}_{\beta }={\Im }^{j+1}\left({z}_{\beta }\right)-{\stackrel{^}{\Im }}_{\beta }^{j+1}$, for β = i, i ± 1 we arrive at

Inserting Equation (29) in Equation (32), we obtain

Using the expressions ${\varpi }_{i-1}={\varpi }_{i}-\ell {\varpi }_{i}^{{}^{\prime }}+\frac{{\ell }^{2}}{2!}{\varpi }^{\left(2\right)}\left({\zeta }_{i}\right)$ and ${\varpi }_{i+1}={\varpi }_{i}+\ell {\varpi }_{i}^{\prime }+\frac{{\ell }^{2}}{2!}{\varpi }^{\left(2\right)}\left({\zeta }_{i}\right)$ in Equation (33), we have

Where ${E}_{i}\left(\ell \right)=\frac{{\ell }^{4}}{3}\left(2{\omega }_{1}-{\omega }_{2}\right){\varpi }_{i}^{j+1}\frac{{d}^{3}{\Im }^{j+1}\left({z}_{i}\right)}{d{z}^{3}}+O\left({\ell }^{6}\right)$. Hence, for the choice of ω1 + ω2 = 1/2, we obtain ${E}_{i}\left(\ell \right)=O\left({\ell }^{4}\right)$.

The matrix representation of Equation (34) is

Where Γ = trid(−σε, 2σε, −σε), $\Lambda =trid\left({\ell }^{2}{\omega }_{1}{r}_{i-1}^{j+1},2{\ell }^{2}{\omega }_{1}{r}_{i}^{j+1},{\ell }^{2}{\omega }_{1}{r}_{i+1}^{j+1}\right)$, $\top ={\left[{e}_{1},{e}_{2},\cdots \phantom{\rule{0.3em}{0ex}},{e}_{N-1}\right]}^{t}$, and $\stackrel{^}{E}={\left[-{E}_{1}\left(\ell \right),-{E}_{2}\left(\ell \right),\cdots \phantom{\rule{0.3em}{0ex}},-{E}_{N-1}\left(\ell \right)\right]}^{t}$.

Following Venkata et al. [30], it can be shown that

Where C is a constant, independent of ℓ and ε.

Theorem 4.4. Let $\Im$(z, t) and ${\stackrel{^}{\Im }}_{i}^{j+1}$ be the solution of Equations (5) and (28) at each grid point (zi, tj+1), respectively. Then, the following uniform error bound holds

$maxi,j|ℑ(zi,tj+1)-ℑ^ij+1|≤C(τ+ℓ2).$

Proof. By combining the result of Lemmas 3.2 and 4.3, the required bound is obtained.

### 4.1. Temporal Richardson extrapolation

In this section, we use the Richardson extrapolation technique to improve the accuracy and order of convergence of the discrete method (28) in the time direction. For that, we consider the tensor product meshes ${\overline{\Upsilon }}^{N,\tau }$ and ${\overline{\Upsilon }}^{N,\tau /2}$, where both the meshes ${\overline{\Upsilon }}^{N,\tau }$ and ${\overline{\Upsilon }}^{N,\tau /2}$ are uniform and identical in spatial direction and uniform in time with step sizes τ and τ/2, respectively. Let ${\overline{\Upsilon }}_{0}^{N,\tau }={\overline{\Upsilon }}^{N,\tau }\cap {\overline{\Upsilon }}^{N,\tau /2}.$ It is clear that ${\overline{\Upsilon }}_{0}^{N,\tau }={\overline{\Upsilon }}^{N,\tau }$. Furthermore, let and be the solutions of the discrete scheme (28). Then we define

Where $\Im$exp is the Richardson extrapolation of $\Im$, which has an enhanced order of convergence by one in time. This technique is known as Richardson extrapolation technique.

Theorem 4.5 (Error After Extrapolation). Assume that NN0 satisfies the assumptions (28). Let $\Im$ be the solution of Equation (5) and $\Im$exp be the extrapolated solution of Equation (28) on the two nested meshes ${\overline{\Upsilon }}^{N,\tau }$ and ${\overline{\Upsilon }}^{2N,\tau /2}$. Then, the new error bound takes the form

$||(ℑ-ℑexp)(zi,tj+1)||≤C((τ)2+ℓ2).$

Proof. We consider the expansion of as

Where $\Im$n is the approximate solution and Γn, for n = 1, 2, is the remainder term. The function ξ is the solution of the problem

To show the convergence of the Richardson technique, we need to formulate the estimate for the remainder term Γn on . It is clear that , where ∂ϒN, τ/n denotes the boundary of ${\overline{\Upsilon }}^{N,\tau /n}.$ In addition, for all $\left({z}_{i},{t}_{j+1}\right)\in {\Upsilon }^{N,\tau /n}$, where n = 1, 2, we have

Then applying the discrete maximum principle, we get $|{\Gamma }^{n}\left({z}_{i},{t}_{j+1}\right)|\le C\left({\left(\tau \right)}^{2}+{\ell }^{2}\right)$.

Then this estimate for Γn, together with the Equations (37) and (38), yields

Equation (40) gives the following error bound for the Richardson extrapolate scheme:

$||(ℑ-ℑexp)(zi,tj+1)||≤C((τ)2+ℓ2).$

## 5. Numerical examples, results, and discussion

Four test examples are presented to illustrate the efficiency of the proposed method. The exact solutions for these problems are unknown for comparison. Therefore, we use the double mesh principle to estimate the error as given in Gupta et al. [31].

$EεN,τ=max(zi,tj+1)∈ϒN,M|(ℑN,τ(zi,tj+1)-ℑ2N,τ/2(zi,tj+1))|,(Before extrapolation),(EεN,τ)extp=max(zi,tj+1)∈ϒN,M|(ℑextpN,τ(zi,tj+1)-ℑextp2N,τ/2(zi,tj+1))|,(After extrapolation),$

Where $\Im$(zi, tj+1) and $\Im$extp(zi, tj+1) denote the numerical solution obtained before and after extrapolation, respectively, in ϒN,M with N mesh intervals in the spatial direction and M mesh intervals in the temporal direction. The rate of convergence is calculated as

$rεN,τ=log2(EεN,τ/Eε2N,τ/2),(Before extrapolation),(rεN,τ)extp=log2((EεN,τ)extp/(Eε2N,τ/2)extp),(After extrapolation).$

For each N and M, the parameter uniform maximum absolute error and uniform rate of convergence are computed using

$EN,τ=maxεEεN,τ,(Before extrapolation),EextpN,τ=maxεEεextpN,τ,(After extrapolation),rN,τ=log2(EN,τt/E2N,τ/2),(Before extrapolation),rextpN,τ=log2(EextpN,τ/Eextp2N,τ/2),(After extrapolation).$

Example 5.1. Consider the following SPBHE:

${∂ℑ∂t-ε∂2ℑ∂z2+ℑ∂ℑ∂z-(1-ℑ)(ℑ-0.5)=0,(z,t)∈ϒ,ℑ(z,0)=z(1-z2),z∈χz,ℑ(0,t)=0=ℑ(1,t),t∈χt.$

Example 5.2. Consider the following SPBHE:

${∂ℑ∂t-ε∂2ℑ∂z2+ℑ∂ℑ∂z=0,(z,t)∈ϒ,ℑ(z,0)=z(1-z2),z∈χz,ℑ(0,t)=0=ℑ(1,t),t∈χt.$

Example 5.3. Consider the following SPBHE:

${∂ℑ∂t-ε∂2ℑ∂z2+ℑ∂ℑ∂z=(1-ℑ)(ℑ-0.5),(z,t)∈ϒ,ℑ(z,0)=sin(πz),z∈χz,ℑ(0,t)=0=ℑ(1,t),t∈χt.$

Example 5.4. Consider the following SPBHE:

${∂ℑ∂t-ε∂2ℑ∂z2+ℑ∂ℑ∂z=0,(z,t)∈ϒ,ℑ(z,0)=sin(πz),z∈χz,ℑ(0,t)=0=ℑ(1,t)=0,t∈χt.$

The computed ${E}_{\epsilon }^{N,\tau },{E}^{N,\tau }$, rN, τ, $\left({E}_{\epsilon }^{N,\tau }\right)extp$, ${E}_{extp}^{N,\tau }$, and ${r}_{extp}^{N,\tau }$ of the considered problems using the proposed scheme are depicted in Tables 14. From these tables, one can observe that the developed scheme converges independent of the perturbation parameter with the order of convergence one before Richardson extrapolation and with the order of convergence two after Richardson extrapolation. In addition, from these Tables, one can observe the comparison of the scheme with the results in the articles [20, 22, 29]. The comparison confirms that the scheme provides a more accurate result than some schemes in Derzie et al. [29], Kaushik and Sharma [20], and Liu et al. [22]. In Figure 1, the numerical solution profile of Examples 5.1 and 5.2 when N = 26, M = 40, and ε = 2−18 are given, which shows that a strong boundary layer is created near z = 1. The effect of ε on the solution profile is shown in Figure 2 and the figure shows as ε → 0 strong boundary layer is created near z = 1. The effect of the time step on the solution profile is depicted in Figure 3. This figure depicts the existence of the boundary layer at z = 1 with time variable t → 1. Figure 4 shows the uniform convergence of the proposed scheme in the error plot for Examples 5.1–5.4.

TABLE 1

Table 1. Comparison of ${E}_{\epsilon }^{N,\tau },{E}^{N,\tau }$, rN, τ, $\left({E}_{\epsilon }^{N,\tau }\right)extp$, ${E}_{extp}^{N,\tau }$, and ${r}_{extp}^{N,\tau }$ for Example 5.1 with results in Liu et al. [22].

TABLE 2

Table 2. Comparison of ${E}_{\epsilon }^{N,\tau },{E}^{N,\tau }$, rN, τ, $\left({E}_{\epsilon }^{N,\tau }\right)extp$, ${E}_{extp}^{N,\tau }$, and ${r}_{extp}^{N,\tau }$ for Example 5.2 with results in Liu et al. [22].

TABLE 3

Table 3. Comparison of ${E}_{\epsilon }^{N,\tau },{E}^{N,\tau }$, rN, τ, $\left({E}_{\epsilon }^{N,\tau }\right)extp$, ${E}_{extp}^{N,\tau }$, and ${r}_{extp}^{N,\tau }$ for Example 5.3 with results in Kaushik and Sharma [20] for τ = 0.001.

TABLE 4

Table 4. Comparison of ${E}_{\epsilon }^{N,\tau },{E}^{N,\tau }$, rN, τ, $\left({E}_{\epsilon }^{N,\tau }\right)extp$, ${E}_{extp}^{N,\tau }$, and ${r}_{extp}^{N,\tau }$ for Example 5.4 with results in Derzie et al. [29].

FIGURE 1

Figure 1. Numerical solution for N = 64, M = 40, and ε = 2−18 (A) Example 5.1 and (B) Example 5.2.

FIGURE 2

Figure 2. Effect of the ε on behavior of the solution with layer formation (A) Example 5.1 and (B) Example 5.2.

FIGURE 3

Figure 3. Numerical solution for N = 64, M = 40, and ε = 2−16 and at different time levels (A) Example 5.1 and (B) Example 5.2.

FIGURE 4

Figure 4. Example 5.1 on (A), Example 5.2 on (B), Example 5.3 on (C), and Example 5.4 on (D), the error plot for different values of ε.

## 6. Conclusion

A higher order parameter uniform numerical scheme for singularly perturbed unsteady Burger-Huxley equation is presented. The developed scheme constitutes the implicit Euler in the time variable and the exponential spline method in the space variable. In addition, Richardson's extrapolation technique is implemented to improve the accuracy in the time direction. The stability and convergence analysis of the scheme are discussed and proved. It is found that the scheme proposed is uniformly convergent with linear and second order convergent before and after Richardson extrapolation. Four test examples are presented to illustrate the efficiency of the proposed scheme. It is found that the proposed scheme gives better results than the methods described in Derzie et al. [29], Kaushik and Sharma [20], and Liu et al. [22].

## 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

ITD and GFD: conceptualization, investigation and formal analysis, and visualization. ITD: software programming and writing—original draft. GFD: writing—review and editing. All authors contributed to the article and approved the submitted version.

## Acknowledgments

The authors would like to thank the editor and referees for careful reading and giving prolific comments.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

1. Duan L, Lu Q. Bursting oscillations near codimension-two bifurcations in the chay neuron model. Int J Nonlinear Sci Numer Simulat. (2006) 7:59–64. doi: 10.1515/IJNSNS.2006.7.1.59

2. Liu S, Fan T, Lu Q. The spike order of the winnerless competition (wlc) model and its application to the inhibition neural system. Int J Nonlinear Sci Num Simulat. (2005) 6:133–38. doi: 10.1515/IJNSNS.2005.6.2.133

3. Aronson DG, Weinberger HF. Multidimensional nonlinear diffusion arising in population genetics. Adv Math. (1978) 30:33–76. doi: 10.1016/0001-8708(78)90130-5

4. Zhang, G.-J., Xu, J.-X., Yao H, Wei, R.-X. Mechanism of bifurcation-dependent coherence resonance of an excitable neuron model. Int J Nonlinear Sci Num Simulat. (2006) 7:447–50. doi: 10.1515/IJNSNS.2006.7.4.447

5. Lewis E, Reiss R, Hamilton H, Harmon L, Hoyle G, Kennedy D. An electronic model of the neuron based on the dynamics of potassium and sodium ion fluxes. Neural Theory and Modeling. Stanford University Press (1964). pp. 154–89.

6. Ismail HN, Raslan K, Abd Rabboh AA. Adomian decomposition method for burger's-huxley and burger's-fisher equations. Appl Math Comput. (2004) 159:291–301. doi: 10.1016/j.amc.2003.10.050

7. Javidi M, Golbabai A. A new domain decomposition algorithm for generalized burger's-huxley equation based on chebyshev polynomials and preconditioning. Chaos Solitons Fractals. (2009) 39:849–57. doi: 10.1016/j.chaos.2007.01.099

8. Estevez P. Non-classical symmetries and the singular manifold method: the burgers and the burgers-huxley equations. J Phys A Math Gen. (1994) 27:2113. doi: 10.1088/0305-4470/27/6/033

9. Krisnangkura M, Chinviriyasit S, Chinviriyasit W. Analytic study of the generalized burger's-huxley equation by hyperbolic tangent method. Appl Math Comput. (2012) 218:10843–7. doi: 10.1016/j.amc.2012.04.044

10. Satsuma J, Ablowitz M, Fuchssteiner B, Kruskal M. Topics in soliton theory and exactly solvable nonlinear equations. Phys Rev Lett. (1987). doi: 10.1142/9789814542210

11. Wang X, Zhu Z, Lu Y. Solitary wave solutions of the generalised burgers-huxley equation. J Phys A Math Gen. (1990) 23:271. doi: 10.1088/0305-4470/23/3/011

12. Hashim I, Noorani MSM, Al-Hadidi MS. Solving the generalized burgers-huxley equation using the adomian decomposition method. Math Comput Model. (2006) 43:1404–11. doi: 10.1016/j.mcm.2005.08.017

13. Hashim I, Noorani M, Batiha B. A note on the adomian decomposition method for the generalized huxley equation. Appl Math Comput. (2006) 181:1439–45. doi: 10.1016/j.amc.2006.03.011

14. Khattak AJ. A computational meshless method for the generalized burger's-huxley equation. Appl Math Model. (2009) 33:3718–29. doi: 10.1016/j.apm.2008.12.010

15. Mohammadi R. B-spline collocation algorithm for numerical solution of the generalized burger's-huxley equation. Num Methods Partial Differ Equ. (2013) 29:1173–91. doi: 10.1002/num.21750

16. Sari M, Gürarslan G. Numerical solutions of the generalized burgers-huxley equation by a differential quadrature method. Math Probl Eng. (2009) 2009:765. doi: 10.1155/2009/370765

17. Appadu AR, Inan B, Tijani YO. Comparative study of some numerical methods for the Burgers-Huxley equation. Symmetry. (2019) 11:1333. doi: 10.3390/sym11111333

18. Appadu AR, Tijani YO, Aderogba AA. On the performance of some NSFD methods for a 2-D generalized Burgers-Huxley equation. J Diff Equ Appl. (2021) 27:1537–73. doi: 10.1080/10236198.2021.1999433

19. Appadu AR, Tijani YO. 1D Generalised burgers-huxley: proposed solutions revisited and numerical solution using FTCS and NSFD methods. Front Appl Math Stat. (2022) 1:773733. doi: 10.3389/fams.2021.773733

20. Kaushik A, Sharma M. A uniformly convergent numerical method on non-uniform mesh for singularly perturbed unsteady burger-huxley equation. Appl Math Comput. (2008) 195:688–706. doi: 10.1016/j.amc.2007.05.067

21. Gupta V, Kadalbajoo MK. A singular perturbation approach to solve burgers-huxley equation via monotone finite difference scheme on layer-adaptive mesh. Commun Nonlinear Sci Num Simulat. (2011) 16:1825–44. doi: 10.1016/j.cnsns.2010.07.020

22. Liu L-B, Liang Y, Zhang J, Bao X. A robust adaptive grid method for singularly perturbed burger-huxley equations. Electron Res Arch. (2020) 28:1439. doi: 10.3934/era.2020076

CrossRef Full Text

23. Kabeto MJ, Duressa GF. Accelerated nonstandard finite difference method for singularly perturbed burger-huxley equations. BMC Res Notes. (2021) 14:1–7. doi: 10.1186/s13104-021-05858-4

24. Jima KM, File DG. Implicit finite difference scheme for singularly perturbed burger-huxley equations. J Partial Diff Equat. (2022) 35:87–100. doi: 10.4208/jpde.v35.n1.6

25. Kehinde OO, Munyakazi JB, Appadu AR. A NSFD discretization of two-dimensional singularly perturbed semilinear convection-diffusion problems. Front Media. (2022) 8:861276. doi: 10.3389/fams.2022.861276

26. Zahra WK, Ashraf MEM. Numerical solution of two-parameter singularly perturbed boundary value problems via exponential spline. J King Saud Univer Sci. (2013) 25:201–8. doi: 10.1016/j.jksus.2013.01.003

27. O'Malley JRE. Introduction to singular perturbations. Technical report, New York Univ Ny Courant Inst of Mathematical Sciences (1974).

28. Nichols N. On the numerical integration of a class of singular perturbation problems. J Optim Theory Appl. (1989) 60:2050034. doi: 10.1007/BF00940347

29. Derzie EB, Munyakazi JB, Gemechu T. A parameter-uniform numericalmethod for singularly perturbed burgers' equation. Comput Appl Math. (2022) 41:1–19. doi: 10.1007/s40314-022-01960-w

30. Venkata RKAS, Palli MMK. A numerical approach for solving singularly perturbed convection delay problems via exponentially fitted spline method. Calcolo. (2017) 54:943–61. doi: 10.1007/s10092-017-0215-6

31. Gupta V, Kadalbajoo MK, Dubey RK. A parameter-uniform higher order finite difference scheme for singularly perturbed time-dependent parabolic problem with two small parameters. Int J Comput Math. (2019) 96:474–99. doi: 10.1080/00207160.2018.1432856

Keywords: singular perturbation problem, burger-huxley equation, implicit Euler method, exponential fitted operator method, Richardson extrapolation

Citation: Daba IT and Duressa GF (2023) Numerical treatment of singularly perturbed unsteady Burger-Huxley equation. Front. Appl. Math. Stat. 8:1061245. doi: 10.3389/fams.2022.1061245

Received: 04 October 2022; Accepted: 18 November 2022;
Published: 04 January 2023.

Edited by:

Vikas Gupta, LNM Institute of Information Technology, India

Reviewed by:

Youssri Hassan Youssri, Cairo University, Egypt
Appanah Rao Appadu, Nelson Mandela University, South Africa
Y. N. Reddy, National Institute of Technology Warangal, India

Copyright © 2023 Daba and Duressa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Imiru Takele Daba, imirutakele@gmail.com