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

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.


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: with the initial and boundary conditions: where Φ x = ða, bÞ, Φ t = ð0, tÞ, and B is the boundary operator. A mathematical model for fðx, t, z, z x Þ = −βz δ z x + γzð1 − z δ Þ in (1) with the above conditions is known as the generalized Burger's-Fisher (gBF) equation and is expressed as follows: and the boundary conditions: 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.
The exact solution of Equation (3) was given by Chen and Zhang [5] as follows: 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 solu-tions, 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- The equation is expressed as follows: The parameters β, γ, and δ are the constants and parameter η ∈ ð0, 1Þ. The initial and boundary conditions are as follows: The exact solution derived by Wang [19], using nonlinear transformations, is reproduced hereunder: where 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.

Advances in Mathematical Physics
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 finitedimensional 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 sixthorder 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 Section 2, firstand 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.

Compact Finite Difference Method
The spatial domain ϕ x = ða, bÞ is divided into uniform mesh with step iteration x i = a + ih, i = 0, 1, 2, ⋯, N, h = ðb − aÞ/N and for time domain ϕ t = ðt 0 , tÞ, with t 0 = 0, a uniform step of size Δt = t j+1 − t j such that t j = t 0 + 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.

Spatial Derivatives of First
Order. The first-order spatial derivatives for CFDM6 at the inner nodes are calculated as follows [42]: 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/3Þð4φ − 1Þ and ψ = ð2/3Þð2 + φÞ. 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!Þh 6 z ð7Þ i : For the value of the derivative at x 0 , x 1 , x N−1 , and x N , one-sided forward and backward schemes have been implemented, which produce following results: The relations (14) and (15) can be represented in the form of a matrix system as 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/3Þ ð1 − τÞ and σ = ð1/3Þð−1 + 10τÞ. This reduces Equation (18) to a tridiagonal system as follows: For the boundary points, one-sided forward and backward schemes have been implemented, which gives the following results: The second-order derivative can be written in the matrix form as

Implementation of CFDM6
By substituting the values of first-order and second-order derivatives in Equations (3) and (8) (ii) Model-II: generalized Burger's-Huxley equation: Advances in Mathematical Physics where L represents the nonlinear differential operator as defined above. In order to solve this system of ODE's from the t j to t j+1 time level, SSP-RK43 is applied using the following operations: By using the initial condition, zðx, tÞ at every required time level can be calculated.

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 problem dz/dt = LðzÞ has a unique solution if LðzÞ satisfies the following conditions: (1) LðzÞ is a real function (2) LðzÞ is well defined and continuous in the domain of t ∈ Φ t and z ∈ ð−∞, ∞Þ (3) There exists a constant called the Lipschitz constant κ such that jLðz, t, ΔtÞ − Lð_ z, t, ΔtÞj ≤ κjz − _ zj, where t ∈ Φ t and z and _ z be any two different points It is clearly seen that LðzÞ 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 every t ∈ Φ t and z, _ z ∈ ð−∞,∞Þ, there exit a constant κ such that Lemma 3. Any single-step method is consistent if ϕðz, t, 0Þ = Lðz, 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 z t ≡ Lz has a unique solution zðtÞ on Φ t and also zðtÞ ∈ C ðp+1Þ Φ t for p ≥ 1. Using Taylor's series expansion about any point t j , where ξ ∈ ðt j , tÞ. Taking t = t j+1 , one gets Thus, the incremental function is defined as It is computed using the approximate value of z j where the exact value zðt j Þ is required. Hence, z j+1 = z j + Δtϕðzðt j Þ, t j , ΔtÞ, j = 0, 1, 2, ⋯, m − 1. To compute the error using Taylor's series, The approximate value using the SSP-RK43 scheme is The following relation is obtained:

Advances in Mathematical Physics
The value of Δtϕðz j , t j , ΔtÞ is obtained from Δtϕðzðt j Þ, t j , ΔtÞ by using the exact approximate value of z j in place of the exact value of zðt j Þ. According to the SSP-RK43, the approximate value of zðt j+1 Þ is obtained as follows: For the above relation, compute the values of zðt j Þ, z ′ ðt j Þ, z′′ðt j Þ ⋯ z p ðt j Þ as follows: Thus, from these computed values taking t = t j , the error term is obtained as follows: Hence, on simplification, In other words, Thus, the given value of p will give the upper bound, and for the computational purpose, the value of L p ðξ j Þ in Equation (37) is replaced with the max | L p ðξ j Þ | in the temporal domain Φ t . The SSP-RK43 as discussed above is rewritten as  Advances in Mathematical Physics The iterated value of z j+1 can be written as Using Taylor's series expansion, the incremental function becomes From the Theorem 1, the proof for convergence is elaborated as follows: 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 L z < κ and L zz < κ 2 /M where M is the upper bound of convergence. For the incremental function, The backward substitution of (38) and its comparison with general Taylor's series [47] gives c 1 = 1/4, c 2 = 1/2, c 3 = 1/4. Hence, these values generate the inequality as

Advances in Mathematical Physics
It is observed that jϕðz j , t j , ΔtÞj satisfies the Lipschitz condition in z j and is a continuous function in Δt. Thus, it is concluded that SSP-RK43 is convergent.

Stability Analysis
The stability analysis of both the models is discussed below by taking nonlinearity coefficient z = m (say), where m = max z, 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.
(2) Model-II: generalized Burger's-Huxley equation: The matrix T 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 whereas on expanding the exponent as a matrix function where I is the identity matrix, For Model-I and Model-II, consider the transformation matrix P such that P −1 TP = D where D is the diagonal matrix; thus, the relation becomes      CFDM6 (Δt = 0:1) 6.4123E-08 6.4126E-08 6.4129E-08 6.4099E-07 6.4102E-07 6.4105E-07 EFD [49] 2.    :

Advances in Mathematical Physics
Taking P −1 z = v in Equation (46), the differential equation becomes Similarly, as discussed above, the solution of Equation (52) is v = exp ðDtÞv 0 , and the recursive relation is In this diagonal matrix, EðDΔtÞ is an approximate matrix of exp ðDΔtÞ. The diagonal elements of the approximated matrix are E j ðη j ΔtÞ. Implementing Equation (25) on the scalar Equation (44), Thus, the method discussed in Equation (25) is absolutely stable if where Re ðηÞ < 0. The stability of the system exclusively depends on the eigenvalues of the coefficient matrix T of the form ∑ 4 m=0 ðTΔtÞ m /m! which should satisfy Equation (54). The necessary conditions that eigenvalues of T should satisfy are given below [47]: Table 6: Comparison of error norms of Example 4 with β = 1, γ = 1, η = 0:001, δ = 2, h = 0:1, and Δt = 0:01.

Method
Error (time) t = 0:05 t = 0:1 t = 1 t = 5 CFDM6 L ∞ 3.2065E-08 6.4129E-08 6.4105E-07 3.1999E-06 L 2 3.0418E-08 6.0835E-08 6.0812E-07 3.0355E-06    (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 Figures 1 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.

Numerical Experiments
The accuracy of compact finite difference scheme is measured using the L 2 and L ∞ error norms, which are defined as follows:     (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],

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 and L 2 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.

Data Availability
The complete data is in the manuscript.

Conflicts of Interest
The authors declare that they have no conflicts of interest.