Front. Appl. Math. Stat. Frontiers in Applied Mathematics and Statistics Front. Appl. Math. Stat. 2297-4687 Frontiers Media S.A. 10.3389/fams.2022.1061245 Applied Mathematics and Statistics Original Research Numerical treatment of singularly perturbed unsteady Burger-Huxley equation Daba Imiru Takele 1 * Duressa Gemechis File 2 1Department of Mathematics, Dilla University, Dilla, Ethiopia 2Department of Mathematics, Jimma University, Jimma, Ethiopia

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

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

This article was submitted to Numerical Analysis and Scientific Computation, a section of the journal Frontiers in Applied Mathematics and Statistics

 04 01 2023 2022 8 1061245 04 10 2022 18 11 2022 Copyright © 2023 Daba and Duressa. 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.

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 O(τ2 + ℓ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.

 singular perturbation problem burger-huxley equation implicit Euler method exponential fitted operator method Richardson extrapolation 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 {£εℑ(z,t)=∂ℑ∂t-ε∂2ℑ∂z2+αℑ∂ℑ∂z-θ(1-ℑ)(ℑ-λ)=0,(z,t)∈ϒ,ℑ(z,0)=ℑ0(z),z∈χ¯z,ℑ(0,t)=Θ0(t),ℑ(1,t)=Θ1(t),t∈χt, Where 0 < ε ≪ 1 is a singular perturbation parameter and α ≥ 1, θ ≥ 0, and λ ∈ (0, 1). The functions Θ0(t), Θ1(t), and ℑ0(z) are assumed to be sufficiently smooth, bounded, and independent of ε.

Lemma 2.1 (Maximum principle). If ℑ∈C2,1(ϒ¯), ℑ|∂ϒ ≥ 0, and £εℑ|ϒ ≥ 0, then ℑ|ϒ¯≥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 ℑ(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 {£εℑ(z,t)=(∂ℑ∂t-ε∂2ℑ∂z2)(z,t)=ζ(z,t,ℑ(z,t),ℑz(z,t)),(z,t)∈ϒ,ℑ(z,0)=ℑ0(z),z∈χ¯z,ℑ(0,t)=Θ0(t),ℑ(1,t)=Θ1(t),t∈(0,T], Where ζ(z,t,ℑ(z,t),ℑz(z,t))=-αℑ∂ℑ∂z+θ(1-ℑ)(ℑ-λ) is the non-linear function of z, t, ℑ(z, t), ℑz(z, t).

To linearize the semi-linear term of Equation (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 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]. ℑ0(z,t)=ℑ0(z)exp(-π2t) Thus, the nonlinear term ζ(z, t, ℑ(z, t), ℑz(z, t)) of Equation (2) can be linearized by applying the Newton-Raphson- Kantorovich approximation approach as ζ(z,t,ℑ(k+1)(z,t),ℑz(k+1)(z,t))≅ζ(z,t,ℑ(k)(z,t),ℑz(k)(z,t))+             (ℑ(k+1)(z,t)-ℑ(k)(z,t))ζℑ|(ℑ,t,ℑ(k)(z,t),ℑz(k)(z,t))+(ℑz(k+1)(z,t)-ℑz(k)(ℑ,t))ζℑz|(z,t,ℑ(k)(z,t),ℑz(k)(z,t))+⋯, Where {ℑ(k)}k=0∞ are a sequence of approximate solutions of ζ(z, t, ℑ(z, t), ℑz(z, t)).

Substituting Equation (4) into Equation (2), we have {£εℑ(k+1)(z,t)=ℑt(k+1)(z,t)-εℑzz(k+1)(z,t)                             +ϖ(z,t)ℑz(k+1)(z,t)+q(z,t)ℑ(k+1)(z,t)                             =v(z,t),ℑ(k+1)(z,0)=ℑ0(z),z∈χ¯z,ℑ(k+1)(0,t)=Θ0(t),t∈χ¯t,ℑ(k+1)(1,t)=Θ1(t),t∈χ¯t, where ϖ(z,t)=-ζuz|(z,t,ℑ(k)(z,t),ℑz(k)(z,t)),q(z,t)=-ζℑ|(z,t,ℑ(k)(z,t),ℑz(k)(z,t)),v(z,t)=ζ(z,t,ℑ(k)(z,t),ℑz(k)(z,t))             -ℑ(k)ζℑ|(z,t,ℑ(k)(z,t),ℑz(k)(z,t))             -ℑz(k)ζℑz|(z,t,ℑ(k)(z,t),ℑz(k)(z,t)).

 3.2. Temporal semi-discretization

We define the uniform mesh for the domain χt as χt=ϒτM={tj=jτ,∀j=0,2,3,⋯,M,τ=1/M} and we approximated the temporal derivative term of Equation (5) using implicit Euler method, which gives {(1+τ£εM)ℑ^k+1(z,tj+1)=-εℑ^zzk+1(z,tj+1)                             +ϖ(z,tj+1)ℑ^zk+1(z,tj+1)                             +q(z,tj+1)ℑ^k+1(z,tj+1)-v(z,tj+1)                             =ℑ^k+1(z,tj),ℑ^k+1(z,0)=ℑ^0k+1(z),z∈χ¯z,ℑ^k+1(0,tj+1)=Θ0(tj+1),0≤j≤M-1,ℑ^k+1(1,tj+1)=Θ1(tj+1),0≤j≤M-1. Since ϖ(z,tj+1)≥ϖ*>0 and q(z,tj+1)≥q*>0,z∈χ¯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 ℑ^k+1(z,tj+1) yields ℑ^k+1(z,tj+1)=ℑ^k+1(z,tj)+τℑ^tk+1(z,tj)+O((τ)2). This implies ℑ^k+1(z,tj+1)-ℑ^k+1(z,tj)τ=ℑ^tk+1(z,tj)+O(τ). Substituting Equation (5) into Equation (7) we have ℑ^k+1(z,tj+1)-ℑ^k+1(z,tj)τ=-(-εℑ^zzk+1(z,tj)+ϖ(z,tj)ℑ^zk+1(z,tj)+q(z,tj)ℑ^k+1(z,tj)-v(z,tj))+O(τ)⇒(1+τ£εM)ℑ^k+1(z,tj+1)-τv(z,tj)+ℑ^k+1(z,tj)=O(τ)2). Subtracting Equation (6) from Equation (8), and the local truncation error LTEj+1=ℑ^k+1(z,tj+1)-ℑ^k+1(z,tj+1) at (j + 1)th is the solution of a BVP (1+τ£εM)LTEj+1=O((τ)2),LTEj+1(0)=0=LTEj+1(1). 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 ||Ej||∞=||∑i=1jLTEi||∞,   j≤T/τ                         ≤||LTE1||∞+||LTE2||∞+||LTE3||∞+⋯+||ej||∞                         ≤c1j(τ)2   (by   Lemma 3.1)                         ≤c1(jτ)(τ)                         ≤c1T(τ) (jτ≤T)                         ≤C(τ), Where c1, C are positive constants independent of ε and τ.

Lemma 3.3. The solution ℑ^k+1(z) of the Equation (6) is bounded by ||∂mℑ^k+1(z,tj+1)∂zm||χ¯z≤C(1+ε-mexp(-(ϖ*(1-z))ε)),                                  0≤m≤4. 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 -εℑ^zzk+1(z,tj+1)+ϖ(z,tj+1)ℑ^zk+1(z,tj+1)+r(z,tj+1)ℑ^k+1(z,tj+1)=φ(z,tj+1),   ℑ^(0,tj+1)=Θ0(tj+1),ℑ^(1,tj+1)=Θ1(tj+1),z∈χz,0≤j≤M-1, Where r(z,tj+1)=(q(z,tj+1)+1τ) and φ(z,tj+1)=v(z,tj+1)+ℑ^k+1(z,tj)τ.

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

The approximation solution ℑ^ij+1 of Equation (10) obtained by the segment S△(z) passing through the points (zi,ℑ^ij+1) and (zi+1,ℑ^i+1j+1). Each mixed spline segment S△j+1(z) has the form [26]. S△j+1(z)=Aieψ(z-zi)+Bie-ψ(z-zi)+Ci(z-zi)+Di,0≤i≤N, 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 S△j+1(zi)=ℑ^ij+1,       S△j+1(zi+1)=ℑ^i+1j+1,d2S△j+1(zi)dz2=Ξi,           d2S△j+1(zi+1)dz2=Ξi+1. Differentiating twice both sides of Equation (11), we get d2S△j+1(z)dz2=ψ2(Aieψ(z-zi)+Bie-ψ(z-zi)). 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 Ai=ℓ2Ξiϱ2-Bi. 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 Ai=e-ϱ(ℓ2Ξi+1ϱ2)-Bie-2ϱ. From Equations (14) and (15), we have Ai=ℓ2(Ξi+1-Ξie-ϱ)2ϱ2sinh(ϱ)  and   Bi=ℓ2(Ξieϱ-Ξi+1)2ϱ2sinh(ϱ). Substituting Equation (12) into Equation (11) at z = zi and z = zi+1, we obtain S△j+1(zi)=ℑ^ij+1=Ai+Bi+Di, and S△j+1(zi+1)=ℑ^i+1j+1                          =Aieϱ+Bie-ϱ+Ciℓ+Di, which implies Ci=(ℑ^i+1j+1-ℑ^ij+1)ℓ-ℓ(Ξi+1-Ξi)ϱ2   andDi=ℑ^ij+1-(ℓ2ϱ2Ξi). The first derivative continuity condition at z = zi, that is S△−1(zi) = S△(zi) gives ℓ2(ω1Ξi-1+2ω2Ξi+ω1Ξi+1)=(ℑ^i-1j+1-2ℑ^ij+1+ℑ^i+1j+1),i=1,2,⋯,N-1, where ω1=(sinh(ϱ)-ϱϱ2sinh(ϱ))  and  ω2=(2ϱcosh(ϱ)-2sinh(ϱ)ϱ2sinh(ϱ)). Now from Equation (10), we have {εΞi-1=ϖj+1(zi-1)ℑ^zj+1(zi-1)+rj+1(zi-1)ℑ^j+1(zi-1)                            -φj+1(zi-1),εΞi=ϖj+1(zi)ℑ^zj+1(zi)+rj+1(zi)ℑ^j+1(zi)-φj+1(zi),εΞi+1=ϖj+1(zi+1)ℑ^zj+1(zi+1)+rj+1(zi+1)ℑ^j+1(zi+1)                            -φj+1(zi+1), Where we approximate ℑ^zj+1(zi-1),ℑ^zj+1(zi) and ℑ^zj+1(zi+1) as {ℑ^zj+1(zi-1)≅-ℑ^i+1j+1+4ℑ^ij+1-3ℑ^i-1j+12ℓ,ℑ^zj+1(zi)≅ℑ^i+1j+1-ℑ^i-1j+12ℓ,ℑ^zj+1(zi+1)≅3ℑ^i+1j+1-4ℑ^ij+1+ℑ^i-1j+12ℓ. Substituting Equation (20) into (19), we have {εΞi-1=ϖj+1(zi-1)(-ℑ^i+1j+1+4ℑ^ij+1-3ℑ^i-1j+12ℓ)                            +rj+1(zi-1)ℑ^j+1(zi-1)-φj+1(zi-1),εΞi=ϖj+1(zi)(ℑ^i+1j+1-ℑ^i-1j+12ℓ)+rj+1(zi)ℑ^j+1(zi)-φj+1(zi),εΞi+1=ϖj+1(zi+1)(3ℑ^i+1j+1-4ℑ^ij+1+ℑ^i-1j+12ℓ)                            +rj+1(zi+1)ℑ^j+1(zi+1)-φj+1(zi+1). Substituting Equation (21) into Equation (18) and rearranging, we get ε(ℑ^i-1j+1-2ℑ^ij+1+ℑ^i+1j+1ℓ2)=(-3ω1ϖi-1j+12ℓ-ω2ϖij+1ℓ+ω1ϖi+1j+12ℓ+ω1ri-1j+1)ℑ^i-1j+1+(2ω1ϖi-1j+1ℓ-2ω1ϖi+1j+1ℓ+2ω2rij+1)ℑ^ij+1+                                                 (-ω1ϖi-1j+12ℓ+ω2ϖij+1ℓ+3ω1ϖi-1j+12ℓ+ω1ri+1j+1)ℑ^i+1j+1=(ω1φi-1j+1+ω2φij+1+ω1φj+1i+1). To grip the effect of the perturbation parameter ε, we multiply the perturbation parameter of Equation (22) by fitting factor σ(ρ), we obtain (-σ(ρ)ρ-3ω1ϖi-1j+12-ω2ϖij+1+ω1ϖi+1j+12+ℓω1ri-1j+1)ℑ^i-1j+1+(2σ(ρ)ρ+2ω1ϖi-1j+1-2ω1ϖi+1j+1+2ℓω2rij+1)ℑ^ij+1+                                                 (-σ(ρ)ρ-ω1ϖi-1j+12+ω2ϖij+1+3ω1ϖi-1j+12+ℓω1ri+1j+1)ℑ^i+1j+1=ℓ(ω1φi-1j+1+ω2φij+1+ω1φi+1j+1), Where ρ=ℓε.

Evaluating the limit of Equation (23) as ℓ → 0 limℓ→0(σ(ρ)ρ)(ℑ^j+1(iℓ-ℓ)-2ℑ^j+1(iℓ)+ℑ^j+1(iℓ+ℓ))+(ω1+ω2)limℓ→0(ϖj+1(iℓ))(ℑ^j+1(iℓ+ℓ)-ℑ^j+1(iℓ-ℓ))=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 ℑ^j+1(z)≈ℑoj+1(z)+ϖj+1(1)ϖj+1(z)(Θ1j+1-ℑ^0j+1(1))exp(-ϖj+1(z)(1-z)ε)+O(ε), Where ℑ0j+1(z) is the solution of the reduced problem ϖj+1(z)∂ℑ^0j+1(z)∂z+rj+1(z)ℑ^0j+1(z)=φj+1(z),                       with ℑ^0j+1(1)=Θ1j+1. Using Taylor's series expansion for ϖj+1(z) about the point z = 1 and restricting to their first terms, Equation (25) becomes ℑ^j+1(z)≈ℑ^0j+1(z)+(Θ1j+1-ℑ^0j+1(1))exp(-ϖj+1(1)(1-z)ε)+O(ε). Equation (26) at zi = iℓ and as ℓ → 0 becomes {limℓ→0ℑ^j+1(iℓ)≈ℑ^0j+1(0)+(Θ1j+1-ℑ^0j+1(1))                             exp(-ϖj+1(1)(1ε-iρ))+O(ε),      limℓ→0ℑ^j+1((i-1)ℓ)≈ℑ^0j+1(0)+(Θ1j+1-ℑ^0j+1(1))                             exp(-ϖj+1(1)(1ε-iρ+ρ))+O(ε),limℓ→0ℑ^j+1((i+1)ℓ)≈ℑ^0j+1(0)+(Θ1j+1-ℑ^0j+1(1))                             exp(-ϖj+1(1)(1ε-iρ-ρ))+O(ε). Plugging the above equations into Equation (24), we get σ(ρ)=ϖj+1(0)ρ(ω1+ω2)coth(ϖj+1(1)ρ2). Finally, from Equations (23) and (27), we get £εN,Mℑ^ij+1=Hij+1,i=1,2,⋯,N-1, 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 |ϑic|≥|ϑi-|+|ϑ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 ℑij+1 satisfies ℑij+1≥0, on ∂ϒ. Then £εN,Mℑij+1≥0 on ϒN,M implies that ℑij+1≥0 at each point of ϒ¯N,M.

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

Lemma 4.2 (Discrete Uniform Stability). The solution ℑ^ij+1 of the discrete problem (28) at (j + 1)th time level and η=min0≤i≤N{rij+1}, where η is some positive constant that satisfies ||ℑ^ij+1||≤||£εN,Mℑ^ij+1||η+max{|ℑ^0j+1|,|ℑ^Nj+1|}. Proof. We define barrier functions (Πij+1)± as (Πij+1)±=Z±ℑ^ij+1, Where Z=||£εN,Mℑ^ij+1||η+max{|ℑ^0j+1|,|ℑ^Nj+1|}.

On the boundary points, we obtain    (Π0j+1)±=Z±ℑ^0j+1=||£εN,Mℑ^ij+1||η+max{|ℑ^0j+1|,|ℑ^Nj+1|}±Θj+1(0)≥0,    (ΠNj+1)±=Z±ℑ^Nj+1=||£εN,Mℑ^ij+1||η+max{|ℑ^0j+1|,|ℑ^Nj+1|}±Θj+1(N)≥0. Now, on the discretized spatial domain ϒℓ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), =±(-σ(ρ)ρ-3ω1ϖi-1j+12-ω2ϖij+1+ω1ϖi+1j+12+ℓω1ri-1j+1)ℑi-1j+1±(2σ(ρ)ρ+2ω1ϖi-1j+1-2ω1ϖi+1j+1+2ℓω2rij+1)ℑ^ij+1±(--σ(ρ)ρ-ω1ϖi-1j+12+ω2ϖij+1+3ω1ϖi-1j+12+ℓω1ri+1j+1)ℑi+1j+1+(ℓω1ri-1j+1+2ℓω2rij+1+ℓω1ri+1j+1)Z                                       ,±ℓ(ω1φi-1j+1+2ω2φij+1+ω1φi+1j+1)+(ℓω1ri-1j+1+2ℓω2rij+1+ℓω1ri+1j+1)Z,    =(ℓω1ri-1j+1+2ℓω2rij+1+ℓω1ri+1j+1)(||£εN,Mℑ^ij+1||η+max{|ℑ^0j+1|,|ℑ^Nj+1|})        ∓ℓ(ω1φi-1j+1+2ω2φij+1+ω1φi+1j+1),                                                         ‘≥0,       since rj+1(zi)≥η>0. By Lemma 4.1, we obtain (Πij+1)±≥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 (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 {ei-1′=dℑ^j+1(zi-1)dz-dℑ^i-1j+1dz=ℓ23d3ℑ^j+1(zi)dz3-ℓ312d4ℑ^j+1(zi)dz4          +ℓ430d5ℑ^j+1(ξi)dz5,ei′=dℑ^j+1(zi)dz-dℑ^ij+1dz=-ℓ26d3ℑ^j+1(zi)dz3-ℓ4120d5ℑ^j+1(ξi)dz5,ei+1′=dℑ^j+1(zi+1)dz-dℑ^i+1j+1dz=ℓ23d3ℑ^j+1(zi)dz3+ℓ312d4ℑ^j+1(zi)dz4          +ℓ430d5ℑ^j+1(ξi)dz5, Where zi−1 < ξ < zi+1.

Substituting σεΞβ=ϖβj+1dℑ^βj+1dz+rβj+1ℑβj+1-φβj+1,  β=i,i ±1 into Equation (18), we get σε(ℑ^i-1j+1-2ℑ^ij+1+ℑ^i+1j+1)          =ℓ2ω1(ϖi-1j+1dℑ^i-1j+1dz+ri-1j+1ℑ^i-1j+1-φi-1j+1)          +2ℓ2ω2(ϖij+1dℑ^ij+1dx+rij+1ℑ^ij+1-φij+1)          +ℓ2ω1(ϖi+1j+1dℑ^i+1j+1dz+ri+1j+1ℑ^i+1j+1-φi+1j+1). Considering the corresponding exact solution to Equation (30), we have σε(ℑ^j+1(zi-1)-2ℑ^j+1(zi)+ℑ^j+1(zi+1))             =ℓ2ω1ϖj+1(zi-1)dℑ^j+1(zi-1)dz+          ℓ2ω1(rj+1(zi-1)ℑ^j+1(zi-1)-φj+1(zi-1))   +2ℓ2ω2(ϖj+1(xi)dℑ^j+1(xi)dz+rj+1(zi)ℑ^j+1(zi))  -2ℓ2ω2gj+1(zi)+ℓ2ω1(ϖj+1(ℓi+1)dℑ^j+1(xi+1)dz+                       rj+1(zi+1)ℑ^j+1(zi+1)-φj+1(zi+1)). Subtracting Equation (30) from Equation (31) and denoting eβ=ℑj+1(zβ)-ℑ^βj+1, for β = i, i ± 1 we arrive at (σε-ℓ2ω1ri-1j+1)ei-1+(-2σε-2ℓ2ω2rij+1)ei                     +(σε-ℓ2ω1ri+1j+1)ei+1=ℓ2(ω1ϖi-1j+1ei-1′+2ω2ϖij+1ei′+ω1ϖi+1j+1ei+1′). Inserting Equation (29) in Equation (32), we obtain (σε-ℓ2ω1ri-1j+1)ei-1+(-2σε-2ℓ2ω2rij+1)ei                     +(σε-ℓ2ω1ri+1j+1)ei+1,   =ℓ43(ω1ϖi-1j+1-ω2ϖij+1+ω1ϖi+1j+1)d3ℑj+1(zi)dz3             +ℓ512(-ω1ϖi-1j+1+ω1ϖi+1j+1)d4ℑj+1(zi)dz4+ℓ660(2ω1ϖi-1j+1-ω2ϖij+1+2ω1ϖi+1j+1)d5ℑj+1(ξi)dz5. Using the expressions ϖi-1=ϖi-ℓϖi ′+ℓ22!ϖ(2)(ζi) and ϖi+1=ϖi+ℓϖi′+ℓ22!ϖ(2)(ζi) in Equation (33), we have    (σε-ℓ2ω1ri-1j+1)ei-1+(-2σε-2ℓ2ω2rij+1)ei+(σε-ℓ2ω1ri+1j+1)ei+1=Ei(ℓ),

Where Ei(ℓ)=ℓ43(2ω1-ω2)ϖij+1d3ℑj+1(zi)dz3+O(ℓ6). Hence, for the choice of ω1 + ω2 = 1/2, we obtain Ei(ℓ)=O(ℓ4).

The matrix representation of Equation (34) is (Γ-Λ)⊤=Ê, Where Γ = trid(−σε, 2σε, −σε), Λ=trid(ℓ2ω1ri-1j+1,2ℓ2ω1rij+1,ℓ2ω1ri+1j+1), ⊤=[e1,e2,⋯,eN-1]t, and Ê=[-E1(ℓ),-E2(ℓ),⋯,-EN-1(ℓ)]t.

Following Venkata et al. [30], it can be shown that ||⊤||≤Cℓ2×O(ℓ4)=C(ℓ2), Where C is a constant, independent of ℓ and ε.

Theorem 4.4. Let ℑ(z, t) and ℑ^ij+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 ϒ¯N,τ and ϒ¯N,τ/2, where both the meshes ϒ¯N,τ and ϒ¯N,τ/2 are uniform and identical in spatial direction and uniform in time with step sizes τ and τ/2, respectively. Let ϒ¯0N,τ=ϒ¯N,τ∩ϒ¯N,τ/2. It is clear that ϒ¯0N,τ=ϒ¯N,τ. Furthermore, let ℑ1(zi,tj+1),for all (zi,tj+1)∈ϒ¯N,τ and ℑ2(zi,tj+1),for all (zi,tj+1)∈ϒ¯N,τ/2 be the solutions of the discrete scheme (28). Then we define ℑexp(zi,tj+1)=2ℑ2(zi,tj+1)-ℑ1(zi,tj+1),(zi,tj+1)∈ϒ¯0N,τ, Where ℑexp is the Richardson extrapolation of ℑ, 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 N ≥ N0 satisfies the assumptions (28). Let ℑ be the solution of Equation (5) and ℑexp be the extrapolated solution of Equation (28) on the two nested meshes ϒ¯N,τ and ϒ¯2N,τ/2. Then, the new error bound takes the form ||(ℑ-ℑexp)(zi,tj+1)||≤C((τ)2+ℓ2). Proof. We consider the expansion of ℑn(zi,tj+1),(zi,tj+1)∈ϒ¯N,τ/n for n=1,2 as ℑn(zi,tj+1)=ℑ(zi,tj+1)-2-(n-1)τξn(zi,tj+1)+Γn(zi,tj+1), Where ℑn is the approximate solution and Γn, for n = 1, 2, is the remainder term. The function ξ is the solution of the problem £εξ(zi,tj+1)=-2-1∂2∂t2ℑ(zi,tj+1),(z,t)∈ϒ                       ξ(z,t)=0,(z,t)∈∂ϒ. To show the convergence of the Richardson technique, we need to formulate the estimate for the remainder term Γn on ϒ¯N,τ/n for n=1,2. It is clear that Γn(zi,tj+1)=0 for all (zi,tj+1)∈∂ϒN,τ/n, where ∂ϒN, τ/n denotes the boundary of ϒ¯N,τ/n. In addition, for all (zi,tj+1)∈ϒN,τ/n, where n = 1, 2, we have |£N,τ/kΓn(zi,tj+1)|=|£N,τ/n(ℑn-ℑ)(zi,tj+1)              -2-(n-1)τ£N,τ/nξ(zi,tj+1)|                            ≤C((τ)2+ℓ2). Then applying the discrete maximum principle, we get |Γn(zi,tj+1)|≤C((τ)2+ℓ2).

Then this estimate for Γn, together with the Equations (37) and (38), yields      (ℑ-ℑexp)(zi,tj+1)=ℑ(zi,tj+1)-(2ℑ2(zi,tj+1)-ℑ1(zi,tj+1),(zi,tj+1))=O((τ)2+ℓ2). 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 ℑ(zi, tj+1) and ℑ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εN,τ,EN,τ, rN, τ, (EεN,τ)extp, EextpN,τ, and rextpN,τ of the considered problems using the proposed scheme are depicted in Tables 1–4. 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.

Comparison of EεN,τ,EN,τ, rN, τ, (EεN,τ)extp, EextpN,τ, and rextpN,τ for Example 5.1 with results in Liu et al. [22].

ε N = 25 26 27 28 29 210
↓ M = 20 40 80 160 320 640
Before richardson extrapolation
2−6 2.8408e-03 1.4682e-03 7.4364e-04 3.7377e-04 1.8731e-04 9.3761e-05
2−8 4.0156e-03 2.0020e-03 9.8557e-04 4.8685e-04 2.4170e-04 1.2040e-04
2−10 4.8909e-03 2.4214e-03 1.1410e-03 5.4081e-04 2.6150e-04 1.2838e-04
2−12 5.3650e-03 2.8201e-03 1.3852e-03 6.4340e-04 2.9416e-04 1.3782e-04
2−14 5.4742e-03 2.9441e-03 1.5118e-03 7.4879e-04 3.5703e-04 1.6336e-04
2−16 5.4827e-03 2.9641e-03 1.5399e-03 7.8120e-04 3.8923e-04 1.9002e-04
2−18 5.4827e-03 2.9644e-03 1.5416e-03 7.8591e-04 3.9638e-04 1.9820e-04
EN, τ 5.4827e-03 2.9644e-03 1.5416e-03 7.8591e-04 3.9638e-04 1.9820e-04
rN, τ 0.8871 0.9433 0.9720 0.9875 0.9999
After richardson extrapolation
2−6 1.1757e-04 3.3841e-05 9.0748e-06 2.3391e-06 6.0263e-07 1.4996e-07
2−8 1.6763e-04 5.5737e-05 1.5609e-05 4.1335e-06 1.0628e-06 2.8261e-07
2−10 2.5551e-04 6.9513e-05 1.8120e-05 4.4961e-06 1.1200e-06 2.8300e-07
2−12 2.7297e-04 7.2989e-05 1.8782e-05 4.5990e-06 1.1281e-06 2.8341e-07
2−14 2.7342e-04 7.3024e-05 1.8789e-05 4.6006e-06 1.1284e-06 2.8346e-07
2−16 2.7377e-04 7.3049e-05 1.8794e-05 4.6016e-06 1.1287e-06 2.8350e-07
2−18 2.7377e-04 7.3049e-05 1.8794e-05 4.6016e-06 1.1287e-06 2.8350e-07
EextpN,τ 2.7377e-04 7.3049e-05 1.8794e-05 4.6016e-06 1.1287e-06 2.8350e-07
rextpN,τ 1.9060 1.9586 2.0301 2.0275 1.9932
Results in Liu et al. [22]
2−6 0.025289 0.017672 0.0090066 0.0048378 0.0025035 0.0012731
2−8 0.038607 0.019497 0.011221 0.0062852 0.0033405 0.0017233
2−10 0.093183 0.070120 0.044773 0.020546 0.010545 0.0051608
2−12 0.17017 0.10083 0.062216 0.039526 0.020493 0.010027
2−14 0.20410 0.16703 0.092110 0.052580 0.028766 0.016523
2−16 0.20450 0.15975 0.12612 0.06.9772 0.038531 0.020526
2−18 0.25614 0.21031 0.13406 0.085618 0.048834 0.026219
EN, τ 0.25614 0.21031 0.13406 0.085618 0.048834 0.026219
rN, τ 0.2844 0.6496 0.6469 0.8100 0.8973

The symbol ↓ shows values of ε in the column.

Comparison of EεN,τ,EN,τ, rN, τ, (EεN,τ)extp, EextpN,τ, and rextpN,τ for Example 5.2 with results in Liu et al. [22].

ε N = 25 26 27 28 29 210
↓ M = 20 40 80 160 320 640
Before richardson extrapolation
2−6 2.2225e-03 1.1579e-03 5.8916e-04 2.9661e-04 1.4877e-04 7.4503e-05
2−8 3.4104e-03 1.7015e-03 8.3475e-04 4.1158e-04 2.0411e-04 1.0162e-04
2−10 4.4223e-03 2.1536e-03 1.0013e-03 4.7090e-04 2.2677e-04 1.1112e-04
2−12 5.0534e-03 2.6284e-03 1.2705e-03 5.7857e-04 2.6065e-04 1.2122e-04
2−14 5.2401e-03 2.8042e-03 1.4305e-03 7.0176e-04 3.2939e-04 1.4772e-04
2−16 5.2884e-03 2.8530e-03 1.4768e-03 7.4538e-04 3.6902e-04 1.7849e-04
2−18 5.3005e-03 2.8654e-03 1.4887e-03 7.5704e-04 3.8034e-04 1.8926e-04
EN, τ 5.3005e-03 2.8654e-03 1.4887e-03 7.5704e-04 3.8034e-04 1.8926e-04
rN, τ 0.8874 0.9447 0.9756 0.9931 1.0069
After richardson extrapolation
2−6 9.4527e-05 2.7189e-05 8.0891e-06 2.3427e-06 6.1922e-07 1.5941e-07
2−8 1.4676e-04 4.7702e-05 1.3418e-05 3.5488e-06 9.9968e-07 2.8086e-07
2−10 2.7265e-04 7.1303e-05 1.7154e-05 4.2520e-06 1.0546e-06 2.6286e-07
2−12 2.7366e-04 7.2313e-05 1.8164e-05 4.3530e-06 1.0975e-06 2.7226e-07
2−14 2.7381e-04 7.2463e-05 1.8371e-05 4.4612e-06 1.1010e-06 2.7378e-07
2−16 2.7384e-04 7.2530e-05 1.8382e-05 4.4819e-06 1.1109e-06 2.7602e-07
2−18 2.7388e-04 7.2533e-05 1.8384e-05 4.4820e-06 1.1110e-06 2.7603e-07
EextpN,τ 2.7388e-04 7.2533e-05 1.8384e-05 4.4820e-06 1.1110e-06 2.7603e-07
rextpN,τ 1.9168 1.9802 2.0362 2.0123 2.0090
Results in Liu et al. [22]
2−6 0.038767 0.018983 0.0096122 0.0050867 0.0026211 0.0013306
2−8 0.044450 0.020109 0.010519 0.059653 0.0031896 0.0016508
2−10 0.083339 0.066120 0.040769 0.019284 0.0091775 0.0049998
2−12 0.18762 0.084106 0.057234 0.03309 0.019041 0.0092260
2−14 0.19755 0.15347 0.083864 0.046945 0.025508 0.014930
2−16 0.21340 0.15016 0.11582 0.061958 0.033441 0.017672
2−18 0.28299 0.17036 0.11805 0.074251 0.041879 0.022413
EN, τ 0.28299 0.17036 0.11805 0.074251 0.041879 0.022413
rN, τ 0.7322 0.5292 0.6689 0.8262 0.9019

The symbol ↓ shows values of ε in the column.

Comparison of EεN,τ,EN,τ, rN, τ, (EεN,τ)extp, EextpN,τ, and rextpN,τ for Example 5.3 with results in Kaushik and Sharma [20] for τ = 0.001.

ε↓N→ 16 32 64 128 256
Before richardson extrapolation
2−6 1.1441e-02 3.3769e-03 1.0928e-03 5.2467e-04 2.5971e-04
2−8 2.4264e-02 1.0286e-02 3.5936e-03 1.2628e-03 5.5544e-04
2−10 2.9305e-02 1.4883e-02 7.0408e-03 3.0236e-03 1.2559e-03
2−12 3.0111e-02 1.6059e-02 8.3119e-03 4.0215e-03 2.1264e-03
2−14 3.0130e-02 1.6165e-02 8.5381e-03 4.3332e-03 1.8652e-03
2−16 3.0130e-02 1.6165e-02 8.5447e-03 4.3660e-03 2.2016e-03
2−18 3.0130e-02 1.6165e-02 8.5447e-03 4.3661e-03 2.2042e-03
2−20 3.0130e-02 1.6165e-02 8.5447e-03 4.3661e-03 2.2042e-03
EN, τ 3.0130e-02 1.6165e-02 8.5447e-03 4.3661e-03 2.2042e-03
rN, τ 0.89833 0.91977 0.96868 0.98609
After richardson extrapolation
2−6 1.1118e-03 4.0166e-04 1.3159e-04 4.1281e-05 1.0687e-05
2−8 1.2086e-03 4.4528e-04 1.6502e-04 4.9481e-05 1.2714e-05
2−10 1.3224e-03 5.3328e-04 1.9675e-04 5.8906e-05 1.5096e-05
2−12 1.4246e-03 5.4330e-04 1.9900e-04 5.9563e-05 1.5227e-05
2−14 1.4275e-03 5.4350e-04 1.9907e-04 5.9565e-05 1.5228e-05
2−16 1.4280e-03 5.4357e-04 1.9909e-04 5.9571e-05 1.5241e-05
2−18 1.4280e-03 5.4357e-04 1.9909e-04 5.9571e-05 1.5241e-05
2−20 1.4280e-03 5.4357e-04 1.9909e-04 5.9571e-05 1.5241e-05
EexpN,τ 1.4280e-03 5.4357e-04 1.9909e-04 5.9571e-05 1.5241e-05
rexpN,τ 1.3935 1.4490 1.7407 1.9667
Results in Kaushik and Sharma [20]
2−6 3.9974e-02 2.1987e-02 1.1611e-02 5.969e-03 3.028e-03
2−8 4.0840e-02 2.2530e-02 1.1904e-02 6.126e-03 3.108e-03
2−10 4.1057e02 2.2675e-02 1.1982e-02 6.166 e-03 3.129e-03
2−12 4.1111e-02 2.2711e-02 1.2002e-02 6.176e-3 3.134e-03
2−14 4.1111e-02 2.2711e-02 1.2002e-02 6.176e-03 3.134e-03
2−16 4.1128e-02 2.2722e-02 1.2008e-02 6.179e-03 3.136e-03
2−18 4.1129e-2 2.2723e-2 1.2008e-2 6.179e-3 3.136e-3
2−20 4.1129e-2 2.2723e-2 1.2008e-2 6.179e-3 3.136e-3
EN, τ 4.1129e-2 2.2723e-2 1.2008e-2 6.179e-3 3.136e-3
rN, τ 0.8560 0.9202 0.9586 0.9784

The symbol ↓ shows values of ε in the column and → shows number nodal points in space direction.

Comparison of EεN,τ,EN,τ, rN, τ, (EεN,τ)extp, EextpN,τ, and rextpN,τ for Example 5.4 with results in Derzie et al. [29].

ε N = 64 128 256 512
↓ M = 20 40 80 160
Before richardson extrapolation
100 1.0681e-02 5.6600e-03 2.9086e-03 1.4738e-03
10−2 1.2760e-02 6.7431e-03 3.4351e-03 1.7144e-03
10−4 1.2850e-02 6.8269e-03 3.5117e-03 1.7800e-03
10−6 1.2850e-02 6.8276e-03 3.5125e-03 1.7808e-03
EN, τ 1.2850e-02 6.8276e-03 3.5125e-03 1.7808e-03
rN, τ 0.91232 0.95888 0.97997
After richardson extrapolation
100 6.5119e-04 1.6178e-04 4.0204e-05 1.0013e-05
10−2 8.1736e-04 1.6222e-04 5.3063e-05 1.3552e-05
10−4 8.9593e-04 2.2944e-04 5.7591e-05 1.4263e-05
10−6 8.9673e-04 2.3022e-04 5.8203e-05 1.4621e-05
EextpN,τ 8.9673e-04 2.3022e-04 5.8203e-05 1.4621e-05
rextpN,τ 1.9617 1.9838 1.9931
Results in Derzie et al. [29]
100 0.01070 0.002336 0.0005212 0.0001290
10−2 0.02550 0.005796 0.001409 0.0003555
10−4 0.09929 0.07841 0.06262 0.04828
10−6 0.09929 0.07841 0.06260 0.04979
EN, τ 0.09929 0.07841 0.06260 0.04979
rN, τ 0.34061 0.32488 0.33031

The symbol ↓ shows values of ε in the column.

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

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

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

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.

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 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. 10.1515/IJNSNS.2006.7.1.59 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. 10.1515/IJNSNS.2005.6.2.133 Aronson DG Weinberger HF. Multidimensional nonlinear diffusion arising in population genetics. Adv Math. (1978) 30:33–76. 10.1016/0001-8708(78)90130-5 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. 10.1515/IJNSNS.2006.7.4.447 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. 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. 10.1016/j.amc.2003.10.050 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. 10.1016/j.chaos.2007.01.099 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. 10.1088/0305-4470/27/6/033 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. 10.1016/j.amc.2012.04.044 Satsuma J Ablowitz M Fuchssteiner B Kruskal M. Topics in soliton theory and exactly solvable nonlinear equations. Phys Rev Lett. (1987). 10.1142/9789814542210 Wang X Zhu Z Lu Y. Solitary wave solutions of the generalised burgers-huxley equation. J Phys A Math Gen. (1990) 23:271. 10.1088/0305-4470/23/3/011 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. 10.1016/j.mcm.2005.08.017 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. 10.1016/j.amc.2006.03.011 Khattak AJ. A computational meshless method for the generalized burger's-huxley equation. Appl Math Model. (2009) 33:3718–29. 10.1016/j.apm.2008.12.010 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. 10.1002/num.21750 Sari M Gürarslan G. Numerical solutions of the generalized burgers-huxley equation by a differential quadrature method. Math Probl Eng. (2009) 2009:765. 10.1155/2009/370765 Appadu AR Inan B Tijani YO. Comparative study of some numerical methods for the Burgers-Huxley equation. Symmetry. (2019) 11:1333. 10.3390/sym11111333 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. 10.1080/10236198.2021.1999433 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. 10.3389/fams.2021.773733 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. 10.1016/j.amc.2007.05.06736353358 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. 10.1016/j.cnsns.2010.07.020 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. 10.3934/era.2020076 Kabeto MJ Duressa GF. Accelerated nonstandard finite difference method for singularly perturbed burger-huxley equations. BMC Res Notes. (2021) 14:1–7. 10.1186/s13104-021-05858-434895302 Jima KM File DG. Implicit finite difference scheme for singularly perturbed burger-huxley equations. J Partial Diff Equat. (2022) 35:87–100. 10.4208/jpde.v35.n1.6 Kehinde OO Munyakazi JB Appadu AR. A NSFD discretization of two-dimensional singularly perturbed semilinear convection-diffusion problems. Front Media. (2022) 8:861276. 10.3389/fams.2022.861276 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. 10.1016/j.jksus.2013.01.003 O'Malley JRE. Introduction to singular perturbations. Technical report, New York Univ Ny Courant Inst of Mathematical Sciences (1974). Nichols N. On the numerical integration of a class of singular perturbation problems. J Optim Theory Appl. (1989) 60:2050034. 10.1007/BF00940347 Derzie EB Munyakazi JB Gemechu T. A parameter-uniform numericalmethod for singularly perturbed burgers' equation. Comput Appl Math. (2022) 41:1–19. 10.1007/s40314-022-01960-w Venkata RKAS Palli MMK. A numerical approach for solving singularly perturbed convection delay problems via exponentially fitted spline method. Calcolo. (2017) 54:943–61. 10.1007/s10092-017-0215-6 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. 10.1080/00207160.2018.1432856