Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 04 January 2023
Sec. Numerical Analysis and Scientific Computation
This article is part of the Research Topic Singularly Perturbed Problems: Asymptotic Analysis and Numerical Solutions View all 6 articles

Numerical treatment of singularly perturbed unsteady Burger-Huxley equation

\nImiru Takele Daba
Imiru Takele Daba1*Gemechis File DuressaGemechis 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

{£ε(z,t)=t-ε2z2+αz-θ(1-)(-λ)=0,(z,t)ϒ,(z,0)=0(z),zχ¯z,(0,t)=Θ0(t),(1,t)=Θ1(t),tχt,    (1)

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-ε2z2)(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],    (2)

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)    (3)

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))+,    (4)

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,    (5)

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),0jM-1,^k+1(1,tj+1)=Θ1(tj+1),0jM-1.    (6)

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(τ).    (7)

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).    (8)

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).    (9)

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(τ),jT/τ.

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

||Ej||=||i=1jLTEi||,   jT/τ                   ||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||χ¯zC(1+ε-mexp(-(ϖ*(1-z))ε)),                                                0m4.

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,0jM-1,    (10)

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 Sj+1(z) has the form [26].

Sj+1(z)=Aieψ(z-zi)+Bie-ψ(z-zi)+Ci(z-zi)+Di,0iN,    (11)

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

Sj+1(zi)=^ij+1,       Sj+1(zi+1)=^i+1j+1,d2Sj+1(zi)dz2=Ξi,           d2Sj+1(zi+1)dz2=Ξi+1.    (12)

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

d2Sj+1(z)dz2=ψ2(Aieψ(z-zi)+Bie-ψ(z-zi)).    (13)

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

d2Sj+1(zi)dz2=Ξi=ψ2(Ai+Bi)

and it yields

Ai=2Ξiϱ2-Bi.    (14)

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

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

d2Sj+1(zi+1)dz2=Ξi+1=ψ2(Aieϱ+Bie-ϱ)

and it gives

Ai=e-ϱ(2Ξi+1ϱ2)-Bie-2ϱ.    (15)

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

Ai=2(Ξi+1-Ξie-ϱ)2ϱ2sinh(ϱ)  and   Bi=2(Ξieϱ-Ξi+1)2ϱ2sinh(ϱ).    (16)

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

Sj+1(zi)=^ij+1=Ai+Bi+Di, and Sj+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).    (17)

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,    (18)

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),    (19)

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.    (20)

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).    (21)

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

ε(^i-1j+1-2^ij+1+^i+1j+12)=(-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).    (22)

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),    (23)

Where ρ=ε.

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

lim0(σ(ρ)ρ)(^j+1(i-)-2^j+1(i)+^j+1(i+))+(ω1+ω2)lim0(ϖj+1(i))(^j+1(i+)-^j+1(i-))=0.    (24)

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(ε),    (25)

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(ε).    (26)

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

{lim0^j+1(i)^0j+1(0)+(Θ1j+1-^0j+1(1))                             exp(-ϖj+1(1)(1ε-iρ))+O(ε),      lim0^j+1((i-1))^0j+1(0)+(Θ1j+1-^0j+1(1))                             exp(-ϖj+1(1)(1ε-iρ+ρ))+O(ε),lim0^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).    (27)

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

£εN,M^ij+1=Hij+1,i=1,2,,N-1,    (28)

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+10, on ∂ϒ. Then £εN,Mij+10 on ϒN,M implies that ij+10 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 η=min0iN{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,0i1. 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||C2,

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,    (29)

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)          +22ω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).    (30)

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))   +22ω2(ϖj+1(xi)d^j+1(xi)dz+rj+1(zi)^j+1(zi))  -22ω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)).    (31)

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σε-22ω2rij+1)ei                     +(σε-2ω1ri+1j+1)ei+1=2(ω1ϖi-1j+1ei-1+2ω2ϖij+1ei+ω1ϖi+1j+1ei+1).    (32)

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

(σε-2ω1ri-1j+1)ei-1+(-2σε-22ω2rij+1)ei                     +(σε-2ω1ri+1j+1)ei+1,   =43(ω1ϖi-1j+1-ω2ϖij+1+ω1ϖi+1j+1)d3j+1(zi)dz3             +512(-ω1ϖi-1j+1+ω1ϖi+1j+1)d4j+1(zi)dz4+660(2ω1ϖi-1j+1-ω2ϖij+1+2ω1ϖi+1j+1)d5j+1(ξi)dz5.    (33)

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σε-22ω2rij+1)ei+(σε-2ω1ri+1j+1)ei+1=Ei(),    (34)

Where Ei()=43(2ω1-ω2)ϖij+1d3j+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

(Γ-Λ)=E^,    (35)

Where Γ = trid(−σε, 2σε, −σε), Λ=trid(2ω1ri-1j+1,22ω1rij+1,2ω1ri+1j+1), =[e1,e2,,eN-1]t, and E^=[-E1(),-E2(),,-EN-1()]t.

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

||||C2×O(4)=C(2),    (36)

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)=22(zi,tj+1)-1(zi,tj+1),(zi,tj+1)ϒ¯0N,τ,    (37)

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 NN0 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),    (38)

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-12t2(zi,tj+1),(z,t)ϒ                               ξ(z,t)=0,(z,t)ϒ.    (39)

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)-(22(zi,tj+1)-1(zi,tj+1),(zi,tj+1))=O((τ)2+2).    (40)

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-ε2z2+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-ε2z2+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-ε2z2+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-ε2z2+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 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
www.frontiersin.org

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

TABLE 2
www.frontiersin.org

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

TABLE 3
www.frontiersin.org

Table 3. 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.

TABLE 4
www.frontiersin.org

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

FIGURE 1
www.frontiersin.org

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

FIGURE 2
www.frontiersin.org

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

FIGURE 3
www.frontiersin.org

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
www.frontiersin.org

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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, yes imirutakele@gmail.com

Disclaimer: 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.