# Numerical Approximation of Generalized Burger’s-Fisher and Generalized Burger’s-Huxley Equation by Compact Finite Difference Method

Authors: Ravneet Kaur; Shallu; Sachin Kumar; V. K. Kukreja
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

In this work, computational analysis of generalized Burger’s-Fisher and generalized Burger’s-Huxley equation is carried out using the sixth-order compact finite difference method. This technique deals with the nonstandard discretization of the spatial derivatives and optimized time integration using the strong stability-preserving Runge-Kutta method. This scheme inculcates four stages and third-order accuracy in the time domain. The stability analysis is discussed using eigenvalues of the coefficient matrix. Several examples are discussed for their approximate solution, and comparisons are made to show the efficiency and accuracy of CFDM6 with the results available in the literature. It is found that the present method is easy to implement with less computational effort and is highly accurate also.


---

## Body

## 1. Introduction

The excerpt approximation of the Navier-Stokes equation is represented by a prominent nonlinear mathematical model known as Burger’s equation. It is the perfect combination of advection and diffusion terms. This equation was introduced by Bateman [1]. Later, Burger [2] extensively worked on this problem, considering the turbulence effect and the statistical aspects. Burger’s equation describes the process of simulating shock wave phenomena, dispersion in a porous medium, heat conduction, diffusion flow, modeling of gas dynamics, traffic flow, propagation and reflection of the nonlinear fluid, boundary layer flow, electrohydrodynamics, sound waves, oil reservoir simulation, etc. The spreading of any species due to the favorable environment of the invasive species or predicting the pattern of spreading was an important issue in the early twenties. The great researcher Fisher [3] proposed a model for the temporal and spatial propagation, depicting the wave of increase in gene frequency in an infinite medium and termed it as Fisher’s equation. It represents the biological processes, ecological systems, pattern formation, etc. Petrovskii and Shigesada [4] combined both the models by assuming that the distribution of species is symmetrical and the environment is homogeneous. The following 1D equation was proposed:
(1)∂z∂t=∂2z∂x2+fx,t,z,zx,inΦ=Φx×Φt,with the initial and boundary conditions:
(2)z=z0,inΦ¯x×t0,Bz=Ω,on∂Φx×Φ¯t,where Φx=a,b, Φt=0,t, and B is the boundary operator. A mathematical model for fx,t,z,zx=−βzδzx+γz1−zδ in (1) with the above conditions is known as the generalized Burger’s-Fisher (gBF) equation and is expressed as follows:
(3)∂z∂t−∂2z∂x2+βzδ∂z∂x−γz1−zδ=0,0≤x≤1,t≥0,subject to the initial condition:
(4)zx,0=12+12tanh−βδx21+δ1/δ,and the boundary conditions:
(5)z0,t=12+12tanhβδ21+δβ2+γ1+δ2β1+δt1/δ,t≥0,(6)z1,t=12+12tanh−βδ21+δ1−β2+γ1+δ2β1+δt1/δ,where β, γ, and δ are the constants. The choice of the value of these constants reduces the model to different forms of PDEs. For γ=0, it reduces to the generalized Burger’s equation. Taking β=0, it becomes the generalized Fisher’s equation. The exact solution of Equation (3) was given by Chen and Zhang [5] as follows:
(7)zx,t=12+12tanh−βδ21+δx−β2+γ1+δ2β1+δt1/δ.Over the past many years, work has been done for the explicit solution of Equation (3). Numerical methods provide a tool for the physical behaviour of the system, although theoretical results are available in the literature. Sari et al. [6] applied the compact finite difference method along with the third-order total variation-diminishing Runge-Kutta scheme in the time domain. Zhao et al. [7] implemented the pseudospectral method using the time discretization by Crank-Nicolson as well as the leapfrog scheme and space discretization by Legendre-Galerkin and Chebyshev-Gauss-Lobatto for nodes. Mohammadi [8] proposed the exponential spline and finite difference approximations. Tatari et al. [9] analyzed the radial basis function collocation technique with the predictor-corrector method. Malik et al. [10] discussed the hybridization of the Exp-function method with the nature-inspired algorithm. Yadav and Jiwari [11] analyzed the finite element analysis with the existence and uniqueness of the weak solution using Galerkin’s finite element method. Macias-Diaz and Gonzalez [12] implemented the finite difference method. Soori [13] obtained the exact solution of the Burger’s-Fisher equation using the variational iteration method and homotopy perturbation method. An exponential time differencing scheme using the method of lines was developed by Bratsos and Khaliq [14]. Gurbuz and Sezer [15] discussed the modified Laguerre matrix-collocation method.The significance and various applications motivated the researchers to compute the analytical and numerical solutions of the Burger’s-Fisher equation. Recently, the dynamical behaviour and exact parametric representations of the traveling wave solutions under different parametric conditions have been discussed by Li [16]. In the findings, the exact monotonic and nonmonotonic kink wave solutions, two-peak solitary wave solutions, and periodic wave solutions, as well as unbounded traveling wave solutions have been obtained. Onyejekwe et al. [17] applied a boundary integral element-based numerical technique, in which the boundary and domain values calculate the fundamental integral inside the domain. The domain integrals due to nonlinearity are considered for computing the solution. Investigation of the global existence and uniqueness of a periodic wave solution has been conducted by Zhang et al. [18].Another important nonlinear equation, describing the interaction between reaction mechanism, convection effect, and diffusion transport is the 1D generalized Burger’s-Huxley (gBH) equation, for whichfx,t,z,zx=−βzδzx+γz1−zδzδ−η. The equation is expressed as follows:
(8)∂z∂t−∂2z∂x2+βzδ∂z∂x=γz1−zδzδ−η,a≤x≤b,t≥0.The parametersβ, γ, and δ are the constants and parameter η∈0,1. The initial and boundary conditions are as follows:
(9)zx,0=η2+η2tanhA1x1/δ,(10)za,t=η2+η2tanhA1a−A2t1/δ,zb,t=η2+η2tanhA1b−A2t1/δ.The exact solution derived by Wang [19], using nonlinear transformations, is reproduced hereunder:
(11)zx,t=η2+η2tanhA1x−A2t1/δ,where
(12)A1=ηδ−β+β2+4γ1+δ4δ+1,A2=βηδ+1−1−η+δ−β+β2+4γδ+12δ+1.Forγ=0, the above model conforms to the generalized Burger’s equation, and considering β=0 and δ=1, the Huxley equation [20] is obtained. For β=0, γ=1, and δ=1, it corresponds to the Fitzhugh-Nagoma equation [21]. Yefimova and Kudryashov [22] applied the Hopf-Cole transformation for solving the gBH equation. The Adomian decomposition method was implemented by Ismail et al. [23]. Gao and Zhao [24] proposed the Exp-function method for a series of exact solutions of the gBH equation. A high-order difference scheme using Taylor’s series expansion was presented by Sari et al. [25]. Celik [26] introduced a numerical method based on the Haar wavelet approach.Zhang et al. [27] reduced the Burger’s-Huxley and Burger’s-Fisher equations into first-order systems and then applied the discontinuous Galerkin method. A numerical scheme based on the finite differences for time integration and cubic B-spline for space integration was proposed by Mohammadi [28]. A fourth-order finite difference method was implemented by Bratsos [29] in a two-time level recurrence relation for the solution of the gBH equation. El-Kady et al. [30] discussed the methods based on cardinal Chebyshev and Legendre basis functions with the Galerkin method, Gauss quadrature formula, and El-Gendi method to convert the problem into ordinary differential equations. Technique based on modified cubic B-spline as the basis function with differential quadrature method was discussed by Singh et al. [31]. The nonstandard finite difference method was analyzed by Zibaei et al. [32]. Bukhari [33] applied local radial basis function differential collocation method. Macias-Diaz [34] used the explicit exponential method. Gilani and Saeed [35] applied the CAS wavelet in conjunction with the Picard technique. Cardinal B-spline wavelet numerical method was used by Shiralashetti and Kumbinarasaiah [36]. A technique based on the hyperbolic-trigonometric tension B-spline method was applied by Alinia and Zarebnia [37]. Loyinmi and Akinfe [38] proposed an algorithm using the coupling of the Elzaki transform with the homotopy perturbation method.Recently, the exact solution has been computed by Kushner and Matviichuk [39] using the theory of finite-dimensional dynamics. Shukla and Kumar [40] applied the numerical scheme based on the Crank-Nicolson finite difference method in collaboration with the Haar wavelet analysis, to obtain the numerical solution. A feed-forward artificial neural network technique is applied by Panghal and Kumar [41] in which the constructed error function is minimized using the quasi-Newton algorithm.Based on the traditional finite difference approximations, Lele [42] proposed well-regulated compact schemes to provide a better representation of shorter proportionate lengths. Many researchers have extended the compact finite difference scheme for linear/nonlinear differential equations, partial differential equations having Dirichlet or Neumann boundary conditions. Ansari et al. [43] implemented the CFD6 scheme for free vibration phenomena of nanobeams in an elastic medium. A similar scheme for incompressible Navier-Stokes and scalar transport equation was analyzed by Boersma [44], a reaction-diffusion equation with delay was approximated by Li et al. [45] and the modified Burger’s equation by Kaur et al. [46].In this work, a numerical scheme based on the sixth-order compact finite difference method (CFDM6) followed by the strong stability-preserving Runge-Kutta method (SSP-RK43) for time integration is used to solve gBF and gBH equations. The advantage of CFDM6 with the SSP-RK43 method is that it computes the results at more mesh points, giving a better approximate solution. The proposed method gives the sixth order of convergence in the spatial domain and the third order in the temporal domain. The proposed method is easy to implement and has less computational cost. The future scope of the method is to solve various arduous linear and nonlinear PDEs.The paper is organized as follows: in Section2, first- and second-order spatial derivatives of the CFDM6 are derived. In Section 3, the proposed method is implemented followed by SSP-RK43. In Section 4, convergence is discussed. In Section 5, stability analysis for the proposed scheme is presented. In Section 6, several test problems are discussed to demonstrate and justify the applicability of the proposed scheme. In Section 7, the conclusion explaining the efficiency of CFDM6 is given.


## 2. Compact Finite Difference Method

The spatial domainϕx=a,b is divided into uniform mesh with step iteration xi=a+ih, i=0,1,2,⋯,N, h=b−a/N and for time domain ϕt=t0,t, with t0=0, a uniform step of size Δt=tj+1−tj such that tj=t0+jΔt, j=0,1,2,⋯, is followed. The method for calculating first-order and second-order derivatives using the compact finite difference scheme is given hereunder.


### 2.1. Spatial Derivatives of First Order

The first-order spatial derivatives for CFDM6 at the inner nodes are calculated as follows [42]:
(13)φz′i−1+z′i+φz′i+1=χzi+2−zi−24h+ψzi+1−zi−12h.For the optimality of the scheme with higher-order accuracy, considerφ=1/3 representing the implicit form of the first-order derivative. The unknown parameters on the other side are calculated by the relation χ=1/34φ−1 and ψ=2/32+φ. By simple calculation, Equation (13) reduces to a sixth-order tridiagonal matrix as a linear system of equations given below with truncation error 4/7!h6zi7:
(14)z′i−1+3z′i+z′i+1=−zi−2−28zi−1+28zi+1+zi+212h,i=2,3,⋯,N−2.For the value of the derivative atx0, x1, xN−1, and xN, one-sided forward and backward schemes have been implemented, which produce following results:
(15)z′0+5z′1=160h−197z0−25z1+300z2−100z3+25z4−3z5,2z′0+11z′1+2z′2=112h−80z0−35z1+136z2−28z3+8z4−z5,2z′N−2+11z′N−1+2z′N=112hzN−5−8zN−4+28zN−3−136zN−2+35zN−1+80zN,5z′N−1+z′N=160h3zN−5−25zN−4+100zN−3−300zN−2+25zN−1+197zN.The relations (14) and (15) can be represented in the form of a matrix system as
(16)Az′=Bz,where
(17)Az′=152112131......131211251z′0z′1z′2..z′N−2z′N−1z′N,Bz=1h−19760−256030060−100602560−360−8012−351213612−2812812−112−112−281202812112−112−281202812112......−112−281202812112−112−281202812112112−8122812−1361235128012360−256010060−30060256019760z0z1z2z3..zN−3zN−2zN−1zN.


### 2.2. Spatial Derivatives of Second Order

Similarly, the second-order derivative is calculated as(18)τz′′i−1+z′′i+τz′′i+1=σzi+2−2zi+zi−24h2+ςzi+1−2zi+zi−1h2.Forτ=0, this equation represents the explicit method to calculate the derivative, and for τ=1/10, it will represent the implicit scheme of the second-order derivative. The unknown constants on the R.H.S. are calculated as ς=4/31−τ and σ=1/3−1+10τ. This reduces Equation (18) to a tridiagonal system as follows:
(19)z′′i−1+10z′′i+z′′i+1=12h2zi+1−2zi+zi−1.For the boundary points, one-sided forward and backward schemes have been implemented, which gives the following results:(20)10z′′0+z′′1=12h211536z0−1555144z1+896z2−77372z3+15136z4−1116z5,z′′N−1+10z′′N=12h211536zN−1555144zN−1+896zN−2−77372zN−3+15136zN−4−1116zN−5.The second-order derivative can be written in the matrix form as(21)ℂz′′=Dz,ℂz′′=10111011101.........11011101110z′′0z′′1z′′2...z′′N−2z′′N−1z′′N,Dz=12h211536−1555144896−7737215136−11161−211−21.........1−211−21−111615136−77372896−155514411536z0z1z2z3.zN−3zN−2zN−1zN.


## 2.1. Spatial Derivatives of First Order

The first-order spatial derivatives for CFDM6 at the inner nodes are calculated as follows [42]:
(13)φz′i−1+z′i+φz′i+1=χzi+2−zi−24h+ψzi+1−zi−12h.For the optimality of the scheme with higher-order accuracy, considerφ=1/3 representing the implicit form of the first-order derivative. The unknown parameters on the other side are calculated by the relation χ=1/34φ−1 and ψ=2/32+φ. By simple calculation, Equation (13) reduces to a sixth-order tridiagonal matrix as a linear system of equations given below with truncation error 4/7!h6zi7:
(14)z′i−1+3z′i+z′i+1=−zi−2−28zi−1+28zi+1+zi+212h,i=2,3,⋯,N−2.For the value of the derivative atx0, x1, xN−1, and xN, one-sided forward and backward schemes have been implemented, which produce following results:
(15)z′0+5z′1=160h−197z0−25z1+300z2−100z3+25z4−3z5,2z′0+11z′1+2z′2=112h−80z0−35z1+136z2−28z3+8z4−z5,2z′N−2+11z′N−1+2z′N=112hzN−5−8zN−4+28zN−3−136zN−2+35zN−1+80zN,5z′N−1+z′N=160h3zN−5−25zN−4+100zN−3−300zN−2+25zN−1+197zN.The relations (14) and (15) can be represented in the form of a matrix system as
(16)Az′=Bz,where
(17)Az′=152112131......131211251z′0z′1z′2..z′N−2z′N−1z′N,Bz=1h−19760−256030060−100602560−360−8012−351213612−2812812−112−112−281202812112−112−281202812112......−112−281202812112−112−281202812112112−8122812−1361235128012360−256010060−30060256019760z0z1z2z3..zN−3zN−2zN−1zN.


## 2.2. Spatial Derivatives of Second Order

Similarly, the second-order derivative is calculated as(18)τz′′i−1+z′′i+τz′′i+1=σzi+2−2zi+zi−24h2+ςzi+1−2zi+zi−1h2.Forτ=0, this equation represents the explicit method to calculate the derivative, and for τ=1/10, it will represent the implicit scheme of the second-order derivative. The unknown constants on the R.H.S. are calculated as ς=4/31−τ and σ=1/3−1+10τ. This reduces Equation (18) to a tridiagonal system as follows:
(19)z′′i−1+10z′′i+z′′i+1=12h2zi+1−2zi+zi−1.For the boundary points, one-sided forward and backward schemes have been implemented, which gives the following results:(20)10z′′0+z′′1=12h211536z0−1555144z1+896z2−77372z3+15136z4−1116z5,z′′N−1+10z′′N=12h211536zN−1555144zN−1+896zN−2−77372zN−3+15136zN−4−1116zN−5.The second-order derivative can be written in the matrix form as(21)ℂz′′=Dz,ℂz′′=10111011101.........11011101110z′′0z′′1z′′2...z′′N−2z′′N−1z′′N,Dz=12h211536−1555144896−7737215136−11161−211−21.........1−211−21−111615136−77372896−155514411536z0z1z2z3.zN−3zN−2zN−1zN.


## 3. Implementation of CFDM6

By substituting the values of first-order and second-order derivatives in Equations (3) and (8), a linear system of equations are obtained for i=0,1,⋯,N:
(i)
Model-I: generalized Burger’s-Fisher equation:(22)∂zi∂t=ℂ−1Dzi−βziδA−1Bzi+γzi1−ziδ≡Lzi.(ii)
Model-II: generalized Burger’s-Huxley equation:(23)∂zi∂t=ℂ−1Dzi−βziδA−1Bzi+γzi1−ziδziδ−η≡Lzi.


### 3.1. SSP-RK43 Scheme

Let(24)dzidt=Lzi,i=0,1,2⋯,N,where L represents the nonlinear differential operator as defined above. In order to solve this system of ODE’s from the tj to tj+1 time level, SSP-RK43 is applied using the following operations:
(25)z1=zj+Δt2Lzj,z2=z1+Δt2Lz1,z3=23zj+13z2+Δt6Lz2,zj+1=z3+Δt2Lz3.By using the initial condition,zx,t at every required time level can be calculated.


## 3.1. SSP-RK43 Scheme

Let(24)dzidt=Lzi,i=0,1,2⋯,N,where L represents the nonlinear differential operator as defined above. In order to solve this system of ODE’s from the tj to tj+1 time level, SSP-RK43 is applied using the following operations:
(25)z1=zj+Δt2Lzj,z2=z1+Δt2Lz1,z3=23zj+13z2+Δt6Lz2,zj+1=z3+Δt2Lz3.By using the initial condition,zx,t at every required time level can be calculated.


## 4. Convergence Analysis

Convergence of the model is investigated below for the desired Equations (22) and (23).Theorem 1.
It is an assumption that the given initial value problemdz/dt=Lz has a unique solution if Lz satisfies the following conditions:
(1)
Lz is a real function(2)
Lz is well defined and continuous in the domain of t∈Φt and z∈−∞,∞(3)
There exists a constant called the Lipschitz constantκ such that Lz,t,Δt−Lż,t,Δt≤κz−ż, where t∈Φt and z and ż be any two different pointsIt is clearly seen thatLz for the generalized Burger’s-Fisher equation and generalized Burger’s-Huxley equation is real, well defined, and continuous. Hence, above theorem is satisfied.Lemma 2.
A single-step method (25) is said to be regular, if the incremental function ϕz,t,Δt satisfies the following conditions:
(1)
The function is well defined and is continuous in the given time and space domain(2)
For everyt∈Φt and z,ż∈−∞,∞, there exit a constant κ such that(26)ϕz,t,Δt−ϕż,t,Δt≤κz−ż.Lemma 3.
Any single-step method is consistent ifϕz,t,0=Lz,t.Theorem 4.
The consistency is the necessary and sufficient condition for the convergence of a regular single-step method with the order (say)p≥1.Proof.
This theorem ensures that the approximate solution converges to the exact solution. For the proof, consider the specific incremental functionϕz,t,Δt. Assume that the given differential equation zt≡Lz has a unique solution zt on Φt and also zt∈Cp+1Φt for p≥1. Using Taylor’s series expansion about any point tj,
(27)zt=ztj+t−tjz′tj+12!t−tj2z′′tj+⋯+1p!t−tjpzptj+1p+1!t−tjp+1zp+1ξj,where ξ∈tj,t. Taking t=tj+1, one gets
(28)ztj+1−ztj=Δtz′tj.
Thus, the incremental function is defined as(29)ϕztj,tj,Δt=Δtz′tj+12!Δt2z′′tj+⋯+1p!Δtpzptj.
It is computed using the approximate value ofzj where the exact value ztj is required. Hence, zj+1=zj+Δtϕztj,tj,Δt,j=0,1,2,⋯,m−1. To compute the error using Taylor’s series,
(30)zj+1=zj+Δtz′j+Δt22!z′′j+Δt33!z′′′j+⋯+Δtpp!zpj+Δtp+1p+1!zp+1ξj.
The approximate value using the SSP-RK43 scheme is(31)zj+1=zj+ΔtLzj+Δt22!L2zj+Δt33!L3zj+⋯+Δtpp!Lpzj.
The following relation is obtained:(32)Δtϕztj,tj,Δt=Δtz′tj+Δt22!L2zj+Δt33!L3zj+⋯+Δtpp!Lpzj.
The value ofΔtϕzj,tj,Δt is obtained from Δtϕztj,tj,Δt by using the exact approximate value of zj in place of the exact value of ztj. According to the SSP-RK43, the approximate value of ztj+1 is obtained as follows:
(33)zj+1=zj+Δtϕzj,tj,Δt+Δt22!ϕ′zj,tj,Δt+Δt33!ϕ′′zj,tj,Δt+⋯.
For the above relation, compute the values ofztj,z′tj,z′′tj⋯zptj as follows:
(34)z′tj=Lztj,tj,z′′tj=Lt+LLz,z′′′tj=Ltt+2′Ltz+L2Lzz+LzLt+LLz.⋮
Thus, from these computed values takingt=tj, the error term is obtained as follows:
(35)Δtp+1p+1!zp+1ξj<ε.
Hence, on simplification,(36)Δtp+1zp+1ξj<εp+1!.
In other words,(37)Δtp+1Lpξj<εp+1!.
Thus, the given value ofp will give the upper bound, and for the computational purpose, the value of Lpξj in Equation (37) is replaced with the max∣Lpξj∣ in the temporal domain Φt. The SSP-RK43 as discussed above is rewritten as
(38)Q1=zj+Δt2Lzj,tj,Q2=Q1+Δt2LQ1,Q3=23zj+13Q2+Δt6LQ2,zj+1=Q3+Δt2LQ3.
The iterated value ofzj+1 can be written as
(39)zj+1=zj+c1Q1+c2Q2+c3Q3.
Using Taylor’s series expansion, the incremental function becomes(40)ϕzj,tj,Δt=Δt−1c1Q1+c2Q2+c3Q3.
From the Theorem1, the proof for convergence is elaborated as follows:
(41)Q1−Q1∗=zj+Δt2Lzj−zj∗+Δt2Lzj∗,Q1−Q1∗≤zj−zj∗+Δt2Lzj−Lzj∗≤1+Δt2κzj−zj∗,Q2−Q2∗=Q1+Δt2LQ1−Q1∗−Δt2LQ1∗,Q2−Q2∗≤Q1−Q1∗+Δt2LQ1−LQ1∗=Q1−Q1∗+Δt2Lzj+Δt2Lzj−Lzj∗+Δt2Lzj∗≤1+Δt2κzj−zj∗+Δt2Lzj+Δt2LzjLzzj+Δt2Lzj2Lzzzj+⋯−Lzj∗−Δt2Lzj∗Lzzj∗−Δt2Lzj∗2Lzzzj∗≤1+Δt2κzj−zj∗+Δt2Lzj−Lzj∗+Δt22LzjLzzj−Lzj∗Lzzj∗+Δt23Lzj2Lzzzj−Lzj∗2Lzzzj∗+⋯≤1+Δtκzj−zj∗+Δt22κ2zj−zj∗=1+Δtκ+Δt2κ2zj−zj∗,Q3−Q3∗=23zj+Q23+Δt2LQ2−23zj∗−Q2∗3−Δt2LQ2∗,Q3−Q3∗=23zj−zj∗+13Q2−Q2∗+Δt2LQ2−LQ2∗≤23zj−zj∗+13Q2−Q2∗+Δt2LQ1+Δt2LQ1−LQ1∗+Δt2LQ1∗≤23zj−zj∗+13Q2−Q2∗+Δt2Lzj−Lzj∗+Δt2LzjLzzj−Lzzj∗≤zj−zj∗+Δt2κ2+Δt2κzj−zj∗+Δt2κzj−zj∗+Δt2κ2zj−zj∗≤zj−zj∗+3Δt2κ+2Δt2κ2zj−zj∗,
As discussed by [47], the free parameters are largely taken according to the range of absolute stability. The other possibility is minimizing the sum of the absolute value of the coefficients of the truncation error. Thus Lz<κ and Lzz<κ2/M where M is the upper bound of convergence. For the incremental function,
(42)ϕzj,tj,Δt−ϕzj∗,tj,Δt=Δt−1c1Q1+c2Q2+c3Q3−c1Q1∗−c2Q2∗−c3Q3∗=Δt−1c1Q1−Q1∗+c2Q2−Q2∗+c3Q3−Q3∗≤Δt−1c11+Δt2κzj−zj∗+c21+Δtκ+Δt2κ2zj−zj∗+c3zj−zj∗+3Δt2κ+2Δt2κ2zj−zj∗≤Δt−1c1+c2+c3+c1+2c2+c3κ2+c2+2c3Δtκ22zj−zj∗.
The backward substitution of (38) and its comparison with general Taylor’s series [47] gives c1=1/4,c2=1/2,c3=1/4. Hence, these values generate the inequality as
(43)ϕzj,tj,Δt−ϕzj∗,tj,Δt≤κ1+12Δtκ+16Δtκ2zj−zj∗.
It is observed thatϕzj,tj,Δt satisfies the Lipschitz condition in zj and is a continuous function in Δt. Thus, it is concluded that SSP-RK43 is convergent.


## 5. Stability Analysis

The stability analysis of both the models is discussed below by taking nonlinearity coefficientz=m (say), where m=maxz, in the entire process to handle the nonlinear term in Equations (22) and (23). The eigenvalue-based technique [45] is followed to establish the stability of the system.
(1)
Model-I: generalized Burger’s-Fisher equation:(44)∂z∂t=ℂ−1Dzi−βmδA−1Bzi+γ1−mδzi,zt=ℂ−1D−βmδA−1B+γ1−mδIzi≡Tzi.(2)
Model-II: generalized Burger’s-Huxley equation:(45)∂zi∂t=ℂ−1Dzi−βmδA−1Bzi+γ1−mδmδ−ηzi,(46)zt=ℂ−1D−βmδA−1Bzi+γ1−mδmδ−ηIzi≡Tzi,dzdt=Tz.The matrixT is constant for both the Model-I and Model-II with the assumption that it has distinct or possibly complex eigenvalues with a negative real part. Using the given initial condition for the analytic solution, the relation becomes
(47)zt=expTtz0,whereas on expanding the exponent as a matrix function where I is the identity matrix,
(48)expTt=I+Tt+Tt22!+Tt33!+⋯.For Model-I and Model-II, consider the transformation matrixP such that P−1TP=D where D is the diagonal matrix; thus, the relation becomes
(49)P−1expTtP=expDt,where
(50)D=η1η2η......ηn−2ηn−1ηn.TakingP−1z=v in Equation (46), the differential equation becomes
(51)dvdt=Dv.Similarly, as discussed above, the solution of Equation (52) is v=expDtv0, and the recursive relation is
(52)vj+1=EDΔtvj.In this diagonal matrix,EDΔt is an approximate matrix of expDΔt. The diagonal elements of the approximated matrix are EjηjΔt. Implementing Equation (25) on the scalar Equation (44),
(53)z′=ηjz.Thus, the method discussed in Equation (25) is absolutely stable if
(54)EjηjΔt<1,where Reη<0. The stability of the system exclusively depends on the eigenvalues of the coefficient matrix T of the form ∑m=04TΔtm/m! which should satisfy Equation (54). The necessary conditions that eigenvalues of T should satisfy are given below [47]:
(i)
For realηj::−2.78<Δtηj<0(ii)
For pure imaginaryηj:−22<Δtηj<22(iii)
For complexηj:Δtηj should lie in the region as given by [48]For different values of parameters, eigenvalues corresponding to gBF and gBH equations are given in Figures1 and 2, respectively. It can be clearly observed that the eigenvalues of all the considered problems satisfy the above defined conditions; therefore, the proposed technique is unconditionally stable.Figure 1
Plot of eigenvalues corresponding to gBF equation withΔt=0.0001 and δ=8.
(a)
β=0.001,γ=0.001(b)
β=1,γ=1(c)
β=0.1,γ=−0.0025Figure 2
Plot of eigenvalues corresponding to gBH equation withΔt=0.0001 and δ=8.
(a)
β=γ=1,η=0.001(b)
β=0.1, γ=0.001, η=0.0001(c)
β=5γ=10, η=0.0001(d)
β=0γ=1, η=0.001


## 6. Numerical Experiments

The accuracy of compact finite difference scheme is measured using theL2 and L∞ error norms, which are defined as follows:
(55)L∞=max0≤i≤Nzi−Zi,L2=h∑i=0Nzi−Zi2,where zi and Zi represent the exact and numerical solutions, respectively, at the node point xi for some fixed time.Example 1.
Consider gBF Equation (3) with the parameters β=0.001 and γ=0.001 for the initial condition as Equation (4) and the boundary conditions as (5) and (6). The exact solution is given by Equation (7). Table 1 gives a comparison of the absolute error for fixed spatial step size h=0.1 and temporal step size Δt=0.0001. Absolute error is calculated at time levels t=0.001,0.010,100 with δ=1 and δ=4. The results are found to be more accurate in comparison to the Adomian decomposition method [23], compact FDM [25], and exponential time differencing method of lines [29]. Figure 3(a) compares numerical and exact solution at different time levels, and Figure 3(b) presents the 3D behaviour of the numerical solution with N=10, Δt=0.01, and δ=8.Example 2.
Consider Equation (3) for β=γ=1 with the initial condition (4) and boundary conditions (5) and (6). The absolute error is compared in Table 2 with those of previous investigators Ismail et al. [23], Sari et al. [25], and Bratsos [29] for h=0.1, Δt=0.0001, and δ=2,8 at t=0.0005 and t=0.0010. Figure 4(a) compares the numerical and exact solution at different time levels, and Figure 4(b) represents the 3D behaviour of numerical solution with N=10, Δt=0.01, and δ=8.Example 3.
Consider Equation (3) for the initial and boundary conditions (4) and (6) with β=0.1 and γ=−0.0025. Table 3 depicts the accuracy of the results obtained by CFDM6, by comparing the absolute error with literature data for h=0.1, Δt=0.0001, and δ=2,4,8. Figure 5(a) compares the numerical and exact solution at different time levels, and Figure 5(b) represents the 3D behaviour of the numerical solution with N=10, Δt=0.01, and δ=8.Example 4.
Consider gBH Equation (8) with the initial and boundary conditions (9) and (10) for parametric values β=γ=1 and η=0.001. The exact solution is given by (11). The absolute error at node points x=0.1,0.5,0.9 is given in Table 4 at t=0.1 and t=1 for h=0.1, Δt=0.0001, and δ=2. Comparison shows that results are better than exponential finite difference scheme [49], hybrid B-spline [50], and tension B-spline collocation method [37]. Table 5 gives a comparison of L∞ error norm for δ=1,4, and 8. Table 6 gives a comparison of L2 and L∞ error norms with δ=2, h=0.1, η=0.001, Δt=0.01 at t=0.05,0.1,1,5. Figure 6(a) represents the absolute error at different time levels with N=10, and Figure 6(b) gives the 3D profile of numerical solution with N=50, Δt=0.01, and δ=8.Example 5.
The gBH Equation (8) is considered for the initial and boundary conditions (9) and (10). The CFDM6 results are evaluated forβ=0.1andγ=0.001, andη=0.0001,h=0.1,Δt=0.0001, andδ=2at timet=0.5andt=0.8are given in Table 7. The absolute error is compared with [37, 49, 50]. The L∞ error norm is compared for CFDM6 with the collocation of cubic B-splines [51], multiscale Runge-Kutta Galerkin method (MGT) [52], and a new kind of tension B-spline function [37] and is presented in Table 8 at different values of δ=1,4,8. Figure 7(a) represents the absolute error at different time levels with N=10, and Figure 7(b) gives the 3D profile of numerical solution with N=50, Δt=0.01, and δ=8.Example 6.
Consider gBH Equation (8) with initial and boundary conditions (9) and (10). The absolute error is compared with the schemes discussed by [37, 49, 50] for β=5, γ=10, η=0.0001, Δt=0.0001, and δ=2 at different node points for time t=0.2,0.5, and 0.8. Tables 9 and 10 give a comparison of absolute error for η=0.0001 and η=0.00001, respectively. Remarkable closeness of numerical and exact solutions can be seen in the tables. Figure 8(a) represents the absolute error at different time levels with N=10, and Figure 8(b) gives the 3D profile of numerical solution with N=50, Δt=0.01, and δ=8.Example 7.
The gBH Equation (8) is subjected to initial and boundary conditions (9) and (10) for β=0, γ=1, and η=0.001. Table 11 compares absolute error of CFDM6 with the Adomian decomposition method (ADM) [23], fourth-order numerical scheme (FDS4) [29], Gauss Chebyshev Galerkin (GCG) [30], and modified cubic B-spline differential quadrature method (MCSDQM) [31] at δ=2, h=0.1, and Δt=0.0001. Table 12 gives the comparison of absolute error for δ=3. The efficiency of the numerical solution to approach the exact solution can be easily seen, and the results are better than those of other methods. Figure 9(a) represents the absolute error at different time levels with N=10, and Figure 9(b) gives the 3D profile of the numerical solution with N=50, Δt=0.01, and δ=8.Table 1
Comparison of absolute error of Example1 with β=0.001, γ=0.001, h=0.1, and Δt=0.0001.
txδ=1δ=4CFDM6ADM [23]CFDM [6]ETDM [14]CFDM6CFDM [6]ETDM [14]0.0010.12.2204E-161.94E-061.01E-071.15E-081.1102E-161.75E-087.71E-090.51.1102E-161.94E-061.04E-073.07E-131.1102E-161.75E-082.07E-130.94.4409E-171.94E-061.01E-071.15E-083.3307E-161.75E-087.71E-090.0100.15.8818E-161.94E-057.53E-076.02E-084.4409E-151.27E-064.05E-080.51.6653E-161.94E-051.04E-068.96E-134.2188E-151.75E-065.56E-130.91.1102E-151.94E-057.53E-076.02E-084.8850E-151.27E-064.05E-081000.12.2204E-16—7.53E-071.01E-075.5511E-16—5.73E-080.51.1102E-15—1.04E-061.50E-112.7756E-15—3.51E-120.91.1102E-16—7.53E-071.01E-071.3323E-15—5.73E-08Figure 3
Graphical representation of solutions corresponding to Example1 with N=10 and Δt=0.001.
(a)
Comparison of numerical and exact solution(b)
Surface plot of numerical solutionTable 2
Comparison of absolute error of Example2 with β=1, γ=1, h=0.1, and Δt=0.0001.
txδ=2δ=8CFDM6ADM [23]CFDM [6]ETDM [14]CFDM6CFDM [6]ETDM [14]0.00050.12.2547E-111.40E-037.62E-055.67E-064.8073E-111.02E-042.44E-060.58.4710E-141.35E-039.14E-055.75E-091.6162E-121.37E-041.82E-100.91.8019E-111.28E-031.02E-045.95E-066.3383E-131.69E-043.15E-060.00100.14.3846E-112.80E-031.50E-041.08E-059.3434E-112.00E-044.65E-060.51.8086E-132.69E-031.83E-041.15E-083.2596E-122.74E-044.02E-100.93.5862E-122.55E-032.00E-041.14E-051.6023E-123.31E-046.00E-06Figure 4
Graphical representation of solutions corresponding to Example2 with N=10 and Δt=0.001.
(a)
Comparison of numerical and exact solution(b)
Surface plot of numerical solutionTable 3
Comparison of absolute error of Example3 with β=0.1, γ=−0.0025, h=0.1, and Δt=0.0001.
txδ=2δ=4δ=8CFDM6CFDM [6]ETDM [14]CFDM6CFDM [6]ETDM [14]CFDM6CFDM [6]ETDM [14]0.10.16.661E-161.21E-059.47E-062.220E-161.34E-056.76E-061.110E-151.47E-054.09E-060.56.661E-162.90E-052.74E-085.551E-163.49E-051.03E-086.661E-163.83E-051.84E-080.92.220E-161.54E-059.57E-067.772E-161.39E-056.92E-083.331E-161.53E-054.24E-060.50.11.341E-161.67E-059.58E-066.661E-162.00E-056.83E-064.441E-162.20E-054.14E-060.51.887E-154.69E-055.18E-082.331E-155.64E-051.93E-086.661E-166.22E-053.47E-080.94.441E-161.71E-059.66E-061.665E-152.07E-057.01E-061.332E-152.28E-054.30E-062.00.15.551E-16—9.59E-061.221E-15—6.86E-061.221E-15—4.20E-060.53.331E-15—5.26E-081.776E-15—1.89E-083.997E-15—3.45E-080.96.661E-16—9.67E-067.772E-16—7.04E-063.331E-16—4.35E-06Figure 5
Graphical representation of solutions corresponding to Example3 with N=10 and Δt=0.001.
(a)
Comparison of numerical and exact solution(b)
Surface plot of numerical solutionTable 4
Comparison of absolute error of Example4 with β=1, γ=1, η=0.001, δ=2, h=0.1, and Δt=0.0001.
Methodt=0.1t=1x=0.1x=0.5x=0.9x=0.1x=0.5x=0.9CFDM6 (Δt=0.1)6.4123E-086.4126E-086.4129E-086.4099E-076.4102E-076.4105E-07EFD [49]2.0510E-065.2339E-062.0511E-063.0562E-068.4901E-063.0564E-06HSCM [50]5.1820E-071.3220E-065.1820E-077.7340E-072.1480E-067.7340E-07UAHT [37]2.8510E-077.8223E-072.8507E-073.0616E-078.5042E-073.0614E-07UAH [37]5.2629E-071.3423E-065.2620E-077.8705E-072.1860E-067.8690E-07UAT [37]5.3131E-071.3585E-065.3121E-077.8706E-072.1861E-067.8691E-07Table 5
Comparison ofL∞ error norm of Example 4 with β=1, γ=1, η=0.001, h=0.1, and Δt=0.001.
Methodt=0.2t=1δ=1δ=4δ=8δ=1δ=4δ=8CFDM6 (Δt=0.1)7.4965E-081.3207E-071.3587E-073.7494E-076.6011E-076.7896E-07MCSCM [51]3.7487E-081.2271E-053.3191E-054.2940E-081.4046E-053.7949E-05MGT [52]4.0305E-081.3193E-053.5687E-054.6849E-081.5325E-054.1407E-05UAHT [37]1.8104E-085.9274E-061.6034E-051.8219E-085.9602E-061.6102E-05UAH [37]4.0069E-081.3118E-053.5485E-054.6833E-081.5321E-054.1400E-05UAT [37]4.0326E-081.3202E-053.5712E-054.6834E-081.5322E-054.1400E-05Table 6
Comparison of error norms of Example4 with β=1, γ=1, η=0.001, δ=2, h=0.1, and Δt=0.01.
MethodError (time)t=0.05t=0.1t=1t=5CFDM6L∞3.2065E-086.4129E-086.4105E-073.1999E-06L23.0418E-086.0835E-086.0812E-073.0355E-06UAHT [37]L∞6.0735E-077.8321E-078.5042E-078.4946E-07L24.4929E-075.7352E-076.2103E-076.2033E-07UAH [37]L∞8.0770E-071.3430E-062.1860E-062.1837E-06L26.2288E-071.0009E-061.5964E-061.5947E-06UAT [37]L∞8.2028E-071.3587E-062.1861E-062.1837E-06L26.3230E-071.0124E-061.5964E-061.5947E-06Figure 6
Error and solution profile of Example4 with N=50 and Δt=0.01.
(a)
Absolute error at different time levels(b)
Surface plot of numerical solutionTable 7
Comparison of absolute error of Example5 with β=0.1, γ=0.001, η=0.0001, δ=2, h=0.1, and Δt=0.0001.
Methodt=0.5t=0.8x=0.1x=0.5x=0.9x=0.1x=0.5x=0.9CFDM6 (Δt=0.1)2.7448E-122.7405E-122.7442E-124.3917E-124.3848E-124.3908E-12EFD [49]4.3493E-111.2069E-104.3494E-114.3758E-111.2154E-104.3759E-11HSCM [50]2.1847E-116.0620E-112.1840E-112.1980E-116.1050E-112.1980E-11UAHT [37]7.3920E-122.0534E-117.3920E-127.3920E-122.0534E-117.3920E-12UAH [37]1.8881E-115.2390E-111.8881E-111.8998E-115.2769E-111.8998E-11UAT [37]1.8892E-115.2427E-111.8892E-111.8999E-115.2772E-111.8999E-11Table 8
Comparison ofL∞ error norm of Example 5 with β=0.1, γ=0.001, η=0.0001, h=0.1, and Δt=0.001 for different values of δ.
Methodt=0.2t=1δ=1δ=4δ=8δ=1δ=4δ=8CFDM6 (Δt=0.1)5.7337E-131.1345E-121.1825E-122.8669E-125.6727E-125.9123E-12MCSCM [51]3.0271E-135.6344E-102.0904E-093.4889E-136.4937E-102.4085E-09MGT [52]3.0804E-135.7325E-102.1267E-093.5806E-136.6634E-102.4720E-09UAHT [37]1.3929E-132.5756E-109.5551E-101.4017E-132.5918E-109.6154E-10UAH [37]3.0631E-135.7006E-102.1148E-093.5847E-136.6629E-102.4718E-09UAT [37]3.0790E-135.7372E-102.1284E-093.5746E-136.6630E-102.4719E-09Figure 7
Error and solution profile of Example5 with N=50 and Δt=0.01.
(a)
Absolute error at different time levels(b)
Surface plot of numerical solutionTable 9
Comparison of absolute error of Example6 with β=5, γ=10, η=0.0001, δ=2, h=0.1, and Δt=0.0001.
txCFDM6 (Δt=0.1)EFD [49]HSCM [50]UAHT [37]0.20.11.2065E-086.58058E-071.971E-078.69755E-080.51.2065E-081.78564E-065.350E-072.41380E-070.91.2065E-086.58087E-071.971E-078.69729E-080.50.13.0158E-087.45354E-072.233E-078.74403E-080.53.0158E-082.06834E-066.198E-072.42887E-070.93.0159E-087.45392E-072.233E-078.74376E-080.80.14.8246E-087.49483E-072.247E-078.74351E-080.54.8247E-082.08190E-066.242E-072.42873E-070.94.8247E-087.49521E-072.247E-078.74324E-08Table 10
Comparison of absolute error of Example6 with β=5, γ=10, η=0.00001, δ=2, h=0.1, and Δt=0.0001.
txCFDM6 (Δt=0.1)EFD [49]HSCM [50]UAHT [37]0.20.11.2066E-102.08154E-086.235E-092.75063E-090.51.2066E-105.64806E-081.692E-087.63381E-090.91.2066E-102.08155E-086.235E-092.75062E-090.50.13.0164E-102.35874E-087.065E-092.76548E-090.53.0164E-106.54514E-081.960E-087.68188E-090.93.0164E-102.35875E-087.065E-092.76547E-090.80.14.8262E-102.37299E-087.108E-092.76547E-090.54.8262E-106.59132E-081.974E-087.68186E-090.94.8262E-102.37300E-087.108E-092.76546E-09Figure 8
Error and solution profile of Example6 with N=50 and Δt=0.01.
(a)
Absolute error at different time levels(b)
Surface plot of numerical solutionTable 11
Comparison of absolute error of Example7 with β=0, γ=1, η=0.001, δ=2, h=0.1, and Δt=0.0001.
txCFDM6ADM [23]FDS4 [29]GCG [30]MCQDQM [31]0.050.13.7491E-085.5890E-071.1176E-064.8110E-074.4924E-070.53.7493E-085.5884E-071.1175E-063.9966E-071.0307E-060.93.7494E-085.5877E-071.1174E-063.9240E-074.4917E-070.100.17.4981E-081.1178E-062.2353E-061.0397E-066.6147E-070.57.4984E-081.1177E-062.2350E-069.5823E-071.7107E-060.97.4987E-081.1175E-062.2347E-069.5091E-076.6139E-071.000.17.4953E-071.1175E-052.2353E-051.1021E-059.9267E-070.57.4956E-071.0074E-052.2350E-051.1057E-052.7793E-060.97.4959E-071.1173E-052.2347E-051.0841E-059.9260E-07Table 12
Comparison of absolute error of Example7 with β=0, γ=1, η=0.001, δ=3, h=0.1, and Δt=0.0001.
txCFDM6ADM [23]FDS4 [29]MCSDQM [31]0.050.13.7499E-081.9841E-063.9673E-061.5946E-060.53.7499E-081.9837E-063.9665E-063.6584E-060.93.7499E-081.9833E-063.9657E-061.5942E-060.100.17.4996E-083.9681E-067.9346E-062.3479E-060.57.4996E-083.9673E-067.9330E-066.0721E-060.97.4996E-083.9665E-067.9314E-062.3475E-061.000.17.4962E-073.9663E-057.9346E-053.5221E-060.57.4962E-073.9655E-057.9330E-059.8610E-060.97.4962E-073.9647E-057.9314E-053.5217E-06Figure 9
Error and solution profile of Example7 with N=50 and Δt=0.01.
(a)
Absolute error at different time levels(b)
Surface plot of numerical solution


## 7. Conclusion

Compact FDM along with the SSP-RK43 scheme has been implemented to solve gBF and gBH equations. Several examples of both the equations are successfully solved with the proposed technique. Absolute error andL2 and L∞ error norms are calculated and compared with the previous results. The results with CFDM6 are found to be better than those with many techniques like the Adomian decomposition method, exponential time differencing method of lines, cubic B-spline collocation method, exponential finite difference scheme, hybrid B-spline collocation, tension B-spline collocation, multiscale Runge-Kutta Galerkin method, and modified cubic B-spline differential quadrature method. Comparison shows that the technique is providing highly accurate results with ease in implementation and less computational effort.



---

*Source: 3346387-2021-11-09.xml*