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.
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: with the initial and boundary conditions: where , , and is the boundary operator. A mathematical model for in (1) with the above conditions is known as the generalized Burger’s-Fisher (gBF) equation and is expressed as follows: subject to the initial condition: 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 , it reduces to the generalized Burger’s equation. Taking , it becomes the generalized Fisher’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 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 which . The equation is expressed as follows:
The parameters , , and are the constants and parameter . The initial and boundary conditions are as follows:
The exact solution derived by Wang [19], using nonlinear transformations, is reproduced hereunder: where
For , the above model conforms to the generalized Burger’s equation, and considering and , the Huxley equation [20] is obtained. For , , and , 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 Section 2, 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 is divided into uniform mesh with step iteration , , and for time domain , with , a uniform step of size such that , , 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]:
For the optimality of the scheme with higher-order accuracy, consider representing the implicit form of the first-order derivative. The unknown parameters on the other side are calculated by the relation and . By simple calculation, Equation (13) reduces to a sixth-order tridiagonal matrix as a linear system of equations given below with truncation error :
For the value of the derivative at , , , and , 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 where
2.2. Spatial Derivatives of Second Order
Similarly, the second-order derivative is calculated as
For , this equation represents the explicit method to calculate the derivative, and for , it will represent the implicit scheme of the second-order derivative. The unknown constants on the R.H.S. are calculated as and . 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
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)Model-I: generalized Burger’s-Fisher equation:(ii)Model-II: generalized Burger’s-Huxley equation:
3.1. SSP-RK43 Scheme
Let where represents the nonlinear differential operator as defined above. In order to solve this system of ODE’s from the to time level, SSP-RK43 is applied using the following operations:
By using the initial condition, 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 problem has a unique solution if satisfies the following conditions: (1) is a real function(2) is well defined and continuous in the domain of and (3)There exists a constant called the Lipschitz constant such that , where and and be any two different points
It is clearly seen that 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 satisfies the following conditions: (1)The function is well defined and is continuous in the given time and space domain(2)For every and , there exit a constant such that
Lemma 3. Any single-step method is consistent if .
Theorem 4. The consistency is the necessary and sufficient condition for the convergence of a regular single-step method with the order (say) .
Proof. This theorem ensures that the approximate solution converges to the exact solution. For the proof, consider the specific incremental function . Assume that the given differential equation has a unique solution on and also for . Using Taylor’s series expansion about any point , where . Taking , one gets Thus, the incremental function is defined as It is computed using the approximate value of where the exact value is required. Hence, . To compute the error using Taylor’s series, The approximate value using the SSP-RK43 scheme is The following relation is obtained: The value of is obtained from by using the exact approximate value of in place of the exact value of . According to the SSP-RK43, the approximate value of is obtained as follows: For the above relation, compute the values of as follows: Thus, from these computed values taking , the error term is obtained as follows: Hence, on simplification, In other words, Thus, the given value of will give the upper bound, and for the computational purpose, the value of in Equation (37) is replaced with the in the temporal domain . The SSP-RK43 as discussed above is rewritten as The iterated value of 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 and where 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 . Hence, these values generate the inequality as It is observed that satisfies the Lipschitz condition in and is a continuous function in . 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 coefficient (say), where , 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:(2)Model-II: generalized Burger’s-Huxley equation:
The matrix 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 is the identity matrix,
For Model-I and Model-II, consider the transformation matrix such that where is the diagonal matrix; thus, the relation becomes where
Taking in Equation (46), the differential equation becomes
Similarly, as discussed above, the solution of Equation (52) is , and the recursive relation is
In this diagonal matrix, is an approximate matrix of . The diagonal elements of the approximated matrix are . Implementing Equation (25) on the scalar Equation (44),
Thus, the method discussed in Equation (25) is absolutely stable if where . The stability of the system exclusively depends on the eigenvalues of the coefficient matrix of the form which should satisfy Equation (54). The necessary conditions that eigenvalues of should satisfy are given below [47]: (i)For real (ii)For pure imaginary (iii)For complex 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.
(a)
(b)
(c)
(a)
(b) , ,
(c) ,
(d) ,
6. Numerical Experiments
The accuracy of compact finite difference scheme is measured using the and error norms, which are defined as follows: where and represent the exact and numerical solutions, respectively, at the node point for some fixed time.
Example 1. Consider gBF Equation (3) with the parameters and 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 and temporal step size . Absolute error is calculated at time levels with and . 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 , , and .
Example 2. Consider Equation (3) for 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 , , and at and . 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 , , and .
Example 3. Consider Equation (3) for the initial and boundary conditions (4) and (6) with and . Table 3 depicts the accuracy of the results obtained by CFDM6, by comparing the absolute error with literature data for , , and . 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 , , and .
Example 4. Consider gBH Equation (8) with the initial and boundary conditions (9) and (10) for parametric values and . The exact solution is given by (11). The absolute error at node points is given in Table 4 at and for , , and . 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 error norm for , and . Table 6 gives a comparison of and error norms with , , , at . Figure 6(a) represents the absolute error at different time levels with , and Figure 6(b) gives the 3D profile of numerical solution with , , and .
Example 5. The gBH Equation (8) is considered for the initial and boundary conditions (9) and (10). The CFDM6 results are evaluated forand, and,,, andat timeandare given in Table 7. The absolute error is compared with [37, 49, 50]. The 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 . Figure 7(a) represents the absolute error at different time levels with , and Figure 7(b) gives the 3D profile of numerical solution with , , and .
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 , , , , and at different node points for time , and . Tables 9 and 10 give a comparison of absolute error for and , 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 , and Figure 8(b) gives the 3D profile of numerical solution with , , and .
Example 7. The gBH Equation (8) is subjected to initial and boundary conditions (9) and (10) for , , and . 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 , , and . Table 12 gives the comparison of absolute error for . 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 , and Figure 9(b) gives the 3D profile of the numerical solution with , , and .
(a) Comparison of numerical and exact solution
(b) Surface plot of numerical solution
(a) Comparison of numerical and exact solution
(b) Surface plot of numerical solution
(a) Comparison of numerical and exact solution
(b) Surface plot of numerical solution
(a) Absolute error at different time levels
(b) Surface plot of numerical solution
(a) Absolute error at different time levels
(b) Surface plot of numerical solution
(a) Absolute error at different time levels
(b) Surface plot of numerical solution
(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 and and 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.
Acknowledgments
Ms. Shallu is thankful to CSIR New Delhi for providing financial assistance in the form of JRF with File No. 09/797(0016)/2018-EMR-I.