A Computational Analysis to Burgers Huxley Equation

Muhammad Saqib,Muhammad Shoaib Arif,Shahid Hasnain and Daoud S.Mashat

1Department of Mathematics,Numl University,Islamabad,44000,Pakistan

2Department of Mathematics,Air University,Islamabad,44000,Pakistan

3Department of Mathematics,Air University,Multan Campus,Multan,66000,Pakistan

4Department of Mathematics,King Abdulaziz University,Jeddah,Saudi Arabia

Abstract:The efficiency of solving computationally partial differential equations can be profoundly highlighted by the creation of precise,higher-order compact numerical scheme that results in truly outstanding accuracy at a given cost.The objective of this article is to develop a highly accurate novel algorithm for two dimensional non-linear Burgers Huxley (BH) equations.The proposed compact numerical scheme is found to be free of superiors approximate oscillations across discontinuities,and in a smooth flow region,it efficiently obtained a high-order accuracy.In particular,two classes of higherorder compact finite difference schemes are taken into account and compared based on their computational economy.The stability and accuracy show that the schemes are unconditionally stable and accurate up to a two-order in time and to six-order in space.Moreover,algorithms and data tables illustrate the scheme efficiency and decisiveness for solving such non-linear coupled system.Efficiency is scaled in terms of L2 and L∞norms,which validate the approximated results with the corresponding analytical solution.The investigation of the stability requirements of the implicit method applied in the algorithm was carried out.Reasonable agreement was constructed under indistinguishable computational conditions.The proposed methods can be implemented for real-world problems,originating in engineering and science.

Keywords:Burgers Huxley equation;finite difference schemes;HOC schemes;Thomas algorithm;Von-Neumann stability analysis

1 Introduction

This paper describes the multiplex schemes solution for two dimensional non-linear Burgers Huxley equation.Such an equation serves as the coupling between theZxx,Zyydiffusive terms andZ(Zx+Zy)the convectional phenomena.This equation is of high importance for showing a prototype model describing the interaction between reaction mechanisms,convection effects and diffusion transports.It is the combination of both Burgers &Huxley phenomena with non-linear term means reactions kind of characteristics behaviour,to capture some features of fluid turbulence which caused by the effects of convection &diffusion [1-3].It is a quantitative paradigm which deals with the flow of electric current through the surface membrane of a giant nerve fibre.Nerve pulse propagation in nerve fibres and wall motion in liquid crystals.Recently research has been measured to investigate two dimensional Burgers Huxley phenomena for understanding the various physical flows in fluid theory [4-6]which leads to implementing a novel methodology for studying new insights [7,8].It is worth mentioning that there is a vast amount of different approaches available in the literature to calculate the solutions of non-linear systems of partial differential equations.Seeking the Burgers Huxley equations numerical solution,wavelet collocation methods for the solution of Burgers Huxley equations [9]have already been studied in combination with variational iteration technique [10,11].Moreover,the propagation of genes (Burger &Fisher) and Reaction-Diffusion (Gray Scott) models [12,13]investigated largely by the technique of computation [14].On the other hand,optimal homotopy asymptotic &homotopy perturbation method was carried out to find the approximate solution of Boussinesq-Burgers equations [15].Finally,some novel techniques also take into account like chaos theory [16],nonlinear optics and fermentation process [17,18].Wazwaz obtained the solitary wave solutions of one dimensional Burgers Huxley equation using tanh-coth method [19].Hashim et al.[20,21]using Adomian Decomposition Method.Molabahrami et al.[22]used the homotopy analysis method to find the solution of one dimensional Burger Huxley equation also Efimova et al.[23]find the travelling wave solution of such equation.Batiha et al.[24]used Hope-Cole transformation with Gao et al.[25]find the exact solution of the generalized Burgers equation.

This research aims to deal with higher-order compact schemes with the finite difference methodology [8].Our primary focus is to attain a compatible scheme which is highly efficient and easy to implement with better accuracy.Although,Burgers Huxley equation can be in three dimensions still some features kept unexplored in the two-dimensional scenario.Let us explorer some new insights in BH equation which consists of the two-dimensional domain which can be written as:

whereZ=Z(l,m,t)is the unknown velocity &(l,m,t)∈Λ×(0T].Laplacian can be defined as

with two dimensional behavior,

alsoPμ,β=(Zμ?1)(β?Z)is a non-linear reaction term.The coefficientξ,ηare advection and reactions coefficients accordingly with 0<β <1 &μ >0.These parameters describe the interaction between reaction mechanisms,convection effects &diffusion transports [26,27].Let us consider the initial condition,

which can be seen from the upcoming Eq.(12).The Dirichlet boundary conditions are given by,

whereΛis a rectangular domain inR2&Z0,p1,p2,q1,q2are given sufficiently smooth functions,andZ(l,m,t)may represent unknown velocity,whereasZ.Zl,Z.Zmrepresents convection terms along with linear diffusionZll,Zmm.Such phenomena perpetuate the ionic mechanisms underlying the initiation and propagation of action potentials in the squid giant axon [28,29].

More generally,it is a challenging task for determining and preservation of physical properties like accuracy,stability,convergence criteria and design efficiency for the given two-dimensional problem.This equation can be an effective procedure for the solution of various deterministic problems in physics,biology and chemical reactions.Also,deals in the investigation of the growth of colonies of bacteria consider population densities or sizes,which are non-negative variables.Most non-linear models of real-life problems are still very challenging to solve either numerically or theoretically.There has recently been much attention devoted to the search for better and more efficient solution methods for determining a solution,analytical or numerical,to non-linear models [30,31].In [31-34]authors present a method used to solve partial equations with the use of artificial neural networks and an adaptive strategy to collocate them.To get the approximate solution of the partial differential equations Deep Neural Networks (DNNs) has been used,which shows impressive results in areas such as visual recognition [35].Recently in [36],authors develop a numerical method with third-order temporal accuracy to solve time-dependent parabolic and first-order hyperbolic partial differential equations.We focused on elaborating further by comparing analytical and numerical techniques.

2 Tanh-Coth Method

The dynamical balance between the non-linear reaction term and diffusive effects which constitute stable waveform after colliding with each other.In (1) the negative coefficients ofZll,ZmmandZ3follow the physical behaviour of two dimensional BH Eq.(1).Such an equation can be converted into the non-linear ordinary differential equation which is as follows:

Letσ=x?et,the wave variable which balances the non-linear reaction term (Pμ,λ(Z))whereμ,λare index values and diffusion transport (the highest derivative involved),we haveM+2=3,MM=1.This enables us to set:PutM=1 in (6) we get

LetY=tanh(γ σ),andσ=((x+y)?et)

Substitutes aforementioned in Eq.(6),we have the following solution to (7)

Arranging the coefficients ofYi,i≥0,and equating these coefficients to zero,the system of algebraic equations ina0,a1,b1,γandeare obtained.By solving the following set of the algebraic system of equations,we have the following form:

In Eq.(10),the solution is of the form:

Case 1:We found thatb1=0,

Case 2:We found thata1=0,

From Cases 1 and 2,the kink solution is of the form:

Now by solving (1) using the tanh-coth method,the analytical solution (kink solution ) is in a compact form in both cases is as follows:

with initial condition:

whereZis the unknown velocity,andγ&σare wavenumbers which are developed during the solution of BH equation.

3 Description of Compact Schemes

Let us discretize the spatial domain which consists ofNandMpositive integers,such thathlandhmpresent step sizes along withlandmdirections,respectively [37].The spatial nodes can be denoted byli,mj,namely,li=ihl,i=0,1,...,N?1,N&mj=jhm,j=0,1,...,M?1,M.For the temporal domain,let us takedtas time-step discretization,τ=T/dt,withtn=Nτ.AlsowithN=0,1,2,...,dt?1 [37-39].Whereτis the temporal step size.SetZZτ={w|w=(w0,w1,w2,...,wdt)T,for anyw∈ZZτ,with some more notations:

for n=0,1,2,...,dt?1.

Implementation Procedure:

Let us we divide (1) into two parts such as:

Now considering one-dimensional steady convection-diffusion equation in the following form:

whereα1,α11are the constants whileβ1,β11are the convective velocities.is the smooth functions oflandmmay represent the reaction,vorticity.Now the three-point scheme is as follows:

Now applying the Taylor series expansion to Eq.(14) we have the following results:

where 0≤n≤dt?1 and the truncation error is

By adding Eqs.(17) and (18) we have the following form (1) which yields:

Apply Crank-Nicholson time discretization,which leads to:

whereT1,T11,T2,T22,ff1,...,are all constants coefficients ofandwhich includesα1,α11,β1,β11,hl,hm,τand constant values.LetZj=[Z1,jZ2,j…ZM?1,j]transposewhere 0≤j≤M.The matrix form of the compact scheme is as follows

By calculating and simplifying the terms,we have the following tridiagonal matrix if of the form:

whereD11matrix is the same as the matrixD33.The Eq.(23) is a tridiagonal block matrix.

The matrix we have generated is diagonally dominant and can solve through Thomas algorithm.Which authenticate the consistency &accuracy of the solution of the formO(h4l+h2m+τ2).

4 Description of Six order Compact Finite Difference Scheme

For complex systems the results will be dependent on the formation of the mesh,We apply higher-order compact scheme at the system in Eq.(1) with a uniform mesh atΔl=Δm=hl=hm.Scheme description is as follows:

Interior Boundary Points:

The compact schemes at interior boundary points are as follows:

Above schemes in Eq.(25) constituteλofamily of tridiagonal structure with parametric valuesForλo=0,we get fourth-order accurate scheme while usingλo=2/11,the scheme becomes sixth-order accurate which leads to[36-39].Also,near boundary points,we have to construct a sixth-order compact scheme to sustain accuracy throughout the two-dimensional domain [36-39].

First Boundary Point 1:

At the first boundary point,the six order compact scheme is of the following form.

Above system in Eq.(26),the coefficients can be found by matching Taylor’s series expansion comparing with various orders up to orderO7,as a result,construction of the linear system is obtained.By constructing the linear system values ofd′s,which can be solved in the usual way to get the following alongldirection,[36,37].Others ones can found in the same way.

2nd Boundary Point 2:

Above system in Eq.(27),by constructing the linear system values ofd′s,which can be solved in the usual way to get the following alongldirection,[36-39].Others ones can found in the same way.

Nth Boundary Point:

At Nth boundary point of six order compact scheme is of the following way:

Above system in Eq.(28),by constructing the linear system values ofd′s,which can be solved in the usual way as done in boundary point 1 and 2.

Implementation Algorithm:

By arranging Eqs.(26)-(28) in the following algorithm:

wherePare mentioned in Eq.(1),also matricesAandBareNm×Nmsparse with triangular nature alongCandDareNl×Nlsparse with triangular in shape.

Theorem:

The truncation error in the compact six order finite difference scheme for equations in the system (1) is,

4.1 Error Analysis

The convergence benchmark,efficiency and accuracy of the proposed scheme in terms of norms can be defined as:

wheredenoted as an analytical solution whilerepresents the numerical solution by mesh points (lo,mq,tn).In this experimentwhereλλis an eigenvalue ofrespectively.

4.2 Stability Analysis

The stability is concerned with the growth or decay of the error produced in the finitedifference solution.For the representation of theoretical analysis,we setP=0 in Eq.(1).Assuming the boundary conditions are accurately propagating,we can apply the Fourier analysis method to our proposed equation.

Definition:For a time-dependent PDE,the corresponding difference scheme is stable in the norm ‖.‖ if there exists a constant M such that

‖en‖≤M‖e0‖,for all n?t≤tF

where M is independent of ?t,?xand initial conditione0.

Following the Von Neumann stability analysis criteria,fix the non-linear terms so that for linear stability,the numerical solution can be displayed in the following way:

whereΓis the amplitude at time leveln,is called the imaginary unit.Φl,Φmleads to wave number inl,andmdirections withΦlhl,Φmhmare phase angles.The amplification factor is defined by

By using Eqs.(30) and (35) and dividing by r.h.s of Eq.(35) and simplifying,we have the following form:

whereR&Sare the compact forms of Eq.(37).For stability,it has to satisfy the following condition:

After simplification to an aforementioned condition which holds true.Therefore,|E|≤1[38-41].Hence the scheme is unconditionally stable.

5 Experimental Results

The novel numerical scheme is compared with the analytical results of Eq.(1) by using tanhcoth method.For this objective,we consider the same parameterα=μ=η=1,and varyingβ.Numerical and analytic solutions are compared and justified in term of error norms to magnify the importance of higher accuracy.

Furthermore,to avoid turbulence,by varyingβvalues in the Tab.1 with grid size(15×15),dt=0.001 and grid space=0.3125 with respect totime=1 is observed.Improvement in accuracy is noted by varying the values ofβparameter.Also,the BH equation produced the best results by using six order compact finite difference scheme.At differentβvalues,Tab.2 indicates error which increased at a very low rate by changing the values ofβfrom high to low which make the comparison to previous work give authentication for accuracy [34].The truncation error is calculated in Tab.3,usingL2,RelativeerrorandL∞with fixed grid size(31×31).By changing time stepsdt=0.001 with the same grid size showed results in the Tab.4.The approximate results using six order compact scheme correspond to error norm are shown in the Tab.6.In this,table the comparison of fourth-order and six order are analyzed by refining the temporal space,which shows this scheme is better than the corresponding fourth-order.In the Tab.7 six order and fourth-order compact finite difference scheme comparison is carried out which measured in term ofL∞norm.Different parameters are also observed under the same scheme.In the Tab.8 scheme efficiency encountered usingL∞,L2&Relativeerrornorms.Graphical representation of numerical schemes on BH equation is observed.Comparison of analytical and numerical results by using fourth and sixth-order compact finite difference scheme has been analyzed.Att=2,β=0.1,dt=0.0001,grid=(21×21)can be seen from the Fig.1.While six order scheme atβ=0.1 with time-spacedt=0.0001 and grid space(21×21)is seen from the Fig.2 which shows more accurate and refine results as compared with Fig.1 using the same parameter.In Figs.3 and 4,analysis shows that the error norm using fourth-order scheme atβ=0.001.While in Fig.5,we chooseβ=0.0001 using a higher-order scheme to analyze error profile at grid size(51×51).In summary,it is aspirant from the figures and tables;the analytical and numerical solutions are best fitted with generation encrypting.In the end,the novel six order compact scheme is the best agreement with the analytical solution.

Comparison between approximation and analytical solutions is made at the final time of computationtime=2sat the critical point(1,1)using fourth-order compact scheme at grid size(15×15).

Table 1:Fourth-order compact scheme

Comparison between approximation and analytical solutions is made at the final time of computationtime=1 at the critical point(1,1)using fourth-order compact scheme at grid size(31×31).

Table 2:Six order compact scheme

Table 3:O(dt2+h4l +h4m)

Table 4:O(dt2+h6l +h6m)

Tab.3 shows error profile data by using fourth-order compact scheme atgridsize=(31×31)for unknown valueZ(l,m,t).Selftime:is the time spent in a function excluding the time spent in its child functions whileTotaltimeis the time to execute the algorithm.

Tab.4 shows error profile data by using Six order compact scheme atgridsize=(31×31)for unknown valueZ(l,m,t).

Tab.5 shows a comparison of two schemes atgridsize=(51×51)andβ=0.001 for unknownsZ(l,m,t).

Tab.6 shows a comparison of two schemes atgridsize=(41×41)for unknownsZ(l,m,t).

Table 5:Error comparison

Table 6:Error comparison

Figure 1:Results obtained by using 4th order scheme at dt=0.00011 &β=0.1 at time=5

Figure 2:Results which obtained by using 6th order scheme at dt=0.0001,&β=0.1 at a time level 2

Figure 3:Results for error estimation by using 4th order scheme at dt=0.001 &β=0.001 at time 1

Figure 4:Results for error estimation by using 6th order scheme at dt=0.001 &β=0.001 at time 1

Figure 5:Results obtained by using 6th order scheme at β=0.0001 at time level 0.5

6 Central Processing Unit Performance

A combinatorial logic circuit executes the mathematical operation for each function in the algorithm within the central processing unit.To establish the platform of CPU performance along physical memory transmission capacity is observed when the higher-order compact scheme is developed by using MATLAB software [35,39,42,43].By increasing the grid size,the number of calculations is increased,and it is difficult to overcome such issue which can take a longer time to execute.Because of numerical schemes efficiency,the computational experiment is done on two different computer machines like Lenovo 6th generation having 2.4 GHz 8 cores and 16 GB memory along 5th generation Dell machine having 4 physical cores and 16 logical cores.Different feathers involved in two computational experiments can be analyzed from the following data tables.

Tab.7 shows results for the different grid using 6th order compact scheme on Lenovo CPU oriented computational machine (MATLAB software).

Tab.8 shows results for the different grid using 6th order compact scheme on DELL CPU oriented computational machine (MATLAB software).

Table 7:Central processing unit performance

Table 8:Central Processing Unit Performance

To check the relative performance and execution time,here in this algorithm which used two machines,Dell and Lenovo,as follows:

where,

The execution time of the Lenovo machine at(17×17)grid size is 170.579s,and Dell machine execution time at the same grid size is 165.123s.To calculate the relative performance we have

Results conclude that Dell machine is 1.033 faster than the Lenovo machine.

Clock cyclecan be defined as:

whereclock cycletime=For CPU time,theclockrate=2GHz,of Lenovo machines and clock rate,is 10s,by increasing theclockratemeans increaseclockcycle.Theclockrateof Lenovo machine=To calculate theclockratewe haveclockcycle=1.2.So

Clock rate performance of Dell machine is,CRDell=

Comparison is performed to analyze Dell with Lenovo machines with both clock rate performance and relative performance.Thus MATLAB handles problems with care,and we can analyze results at each point of the loop and any iteration during computations.

7 Conclusion

Higher-order schemes for determining the two dimensional Burgers Huxley equation was developed in this paper.As it was not studied before by using such schemes of diffusive dissipation of errors.We came to know that the BH equation in two dimensional which is studied to find efficiency,accuracy and stability and by comparing with analytical and numerical approaches in terms ofL2,L∞&relative errors.It is evident from the fact that computed numerical experiments of two dimensional Burgers Huxley equation,solutions obtained by fourth and six order schemes are in good agreement with the analytical solutions.Figures and tables clearly show the tendency of fast and monotonic convergence of the results toward the analytical solution.Also,the computational discretization of the proposed model results in a sparse tridiagonal structure of the matrix,which can be overcome by the Thomas algorithm.Results lead to a remarkable improvement in accuracy,efficiency and computer performance which can be seen from data tables.

Acknowledgement:The authors are thankful to the anonymous referee for their suggestions and helpful comments that improved this article.We are also grateful to Vice-Chancellor,Air University,Islamabad,for providing an excellent research environment and facilities.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.