 American Journal of Computational Mathematics, 2011, 1, 152-158 doi:10.4236/ajcm.2011.13017 Published Online September 2011 (http://www.SciRP.org/journal/ajcm) Copyright © 2011 SciRes. AJCM A Fourth Order Improved Numerical Scheme for the Generalized Burgers—Huxley Equation Athanassios G. Bratsos Department of Mathematics, Technological Educational Institution of Athens, Athens, Greece E-mail: bratsos@teiath.gr Received March 30, 2011; revised May 5, 2011; accepted May 22, 2011 Abstract A fourth order finite-difference scheme in a two-time level recurrence relation is proposed for the numerical solution of the generalized Burgers-Huxley equation. The resulting nonlinear system, which is analyzed for stability, is solved using an improved predictor-corrector method. The efficiency of the proposed method is tested to the kink wave using both appropriate boundary values and conditions. The results arising from the experiments are compared with the relevant ones known in the available bibliography. Keywords: Burgers-Huxley; Finite-Difference Method; Modified Predictor-Corrector 1. Introduction A. Hodgkin and A. Huxley  proposed a model, known henceforth as the Huxley equation, in order to explain the ionic mechanisms underlying the initiation and pro- pagation of action potentials in the squid giant axon. The most general form of the Huxley equation, known as the generalized Burgers-Huxley equation (BgH) [2,3], has the form   =1;0 1,>0txxxuuuu uuuxt  , (1.1) where is a sufficiently often differentiable function, =,uuxt a real parameter, 0, and 0,1>0. Equation (1.1), which models the interaction be-tween reaction mechanisms, convection effects and dif-fusion transport, is the modified Burgers equation for =0 (see  and the references therein), is also the Huxley equation  for =0, =1 and is the Fitz-hugh-Nagoma equation  for =0. Many researchers have used various methods to solve the BgH equation. A theoretical study of the BgH equa-tion was found in Wang et al. , while analytical solu-tions using various techniques in [7-11], etc., have been proposed. As far as the numerical methods are concerned among others the Adomian decomposition method was used by Ismail et al.  for the BgH and the Bur-gers-Fisher equation, and by Hashim et al.  for the BgH equation. Javidi  used the pseudospectral me- thod, while Javidi , Javidi and Golbabai  the spectral collocation method. Batiha et al.  used the variational iteration method and Khattak  the collo-cation method with radial basis functions. Babolian and Saeidian  used the homotopy analysis method, etc. The initial condition associated with Equation (1.1) will be ,0=; 01.uxf xx (1.2) Theoretical Solution It is known  that Equation (1.1) has the following kink wave solution  1,= tanh22uxtkx ct (1.3) in which 241=41k   and  214=121c   1 are the wave number and the velocity respectively. 2. The numerical Method 2.1. Development of the Method 2.1.1. Grid and Solution Vecto r To obtain numerical solutions the region =,Rxt 0< <10,xT with its boundary consisting of R A. G. BRATSOS 153the lines , and is covered with a rectangular mesh of points, , with co-ordinates =0x=1x=0tG,xt =,mnxtn=,mh,mn with . The theoretical solution of Equation (1.1) at the typical mesh point =0,1, ,1mNxtUt will be denoted by and the relevant of an approximating difference scheme by . nmunmU=ntt1, .nNULet the solution vector at time level be 01,,nU==UnnU (2.1) 2.1.2. Boundaries The following were used: 1) The space derivatives at the left boundary were replaced with second order finite-difference replace- ments of the form ( p. 17) =0x2nUU2has014nn=0xx1=32uUh0,h (2.2) 24nUnd 1,u0ash23nU 01=tgt015nn=g21=2uUh0,ut=0xx xhUa (2.3) and with analogous replacements to the right boundary . =1x2) The boundary conditions t (2.4) =0,1=0;un>xx0,t=1 (2.5) were used, while at the other interior points of the grid G the well-known approximants based on the central-dif- ference formulas. 2.2. The Proposed Method Applying Equation (1.1) at each point of the grid G at time level ; leads to a first-order initial-value problem, which is written in a matrix-vector form as ==ntt ,2,n 1 ,Nttfx U0=0 , =Dt t UU=AfxUU01,BtU=fx;>0f(2.6) in which =diagd dDt, == =d gnnt =diagm,nmUia (2.7) gnnmnnmmUU;== =d =diag1tia    (2.8) for are diagonal matrices, =0,1, ,1mN23 4110 1...1=,1 01214 3254112 1.. .1=12114 52AhBh (2.9) or 22 210 1.. .1=,1 012222212 1...1=12 122AhBh (2.10) tridiagonal matrices arising from the use of the boundary values (2.2) - (2.3) or the boundary condition (2.5) re-spectively and f the vector of the initial condition, all of order 2N. Relation (2.6) gives =DAB (2.11) hence can easily be obtained. 2DUsing the recurrence relation =exp ; =0,,,tDttUU  (2.12) where DtU is given by (2.6) and replacing the ma-trix-exponential term with the fourth order rational ap-proximant ( p. 134) gives 222211 21211=.212IDD tID DtUU (2.13) Equation (2.13) using the notations (2.7) - (2.8) and Equation (2.11) leads to the following nonlinear system Copyright © 2011 SciRes. AJCM 154 A. G. BRATSOS  112221 21111 1111 12222121 12 1=21 12 nnnnnn nnnn nnnnnnnn nntABtAB ABABABABttABtAB ABABAB       UUU U .nn nABt UnU(2.14) Let 1=4rh, 22=2rh, 3=2r, 22 24=48rh, 245=12rh, 226=1r2, 237=24rh, 28=2rh4 and 229=12rh. Equation (2.14), when applied to the general mesh point of the grid G, gives 1111 11111211111 11 134121211111 11152 1211 11126 7 2 46 4nnnnnnnmmmm mmmnnnn nn nmmmm mm mnnnnn nmmmmm mnn nnmm mmUrUU rU UUrU rUUUrU UUUU rUr           112111 1111128 1111111 11117121 11111 11112119 1111 22 2 2nnmmnnn nnnmmm mmmnnnnn nnmmmmm mmnnnn nnnmmmm mmmnnmmUUUU rUUr UUUUUrUU       1111111181 1111 191 1 2nnnnnnmmmm m mnnn nmmm mUr UUrUUU    11 11121 134121211252112672112811=2 464 22 nnn nnnnmmmmm mmnnn nnnnnnmmm mmmmmmnnnnnn nmmmmmm mnnnn nmmmm mnnnmmmUrU UrUUUrU rUUUrUUUUUrUrUUUUrU        117 1211 12119111181 191 1 2 2 2; nn nnmm mmnnnnnnnnmmmm mmmmnnnnnnnnnnmmmm mmmmmmnnn nmmm mUr UUU UUrUU UrUUrU UU        = 0,1,...,1.mNnm (2.15) Stability Analysis Following the Fourier method of analysing stability ( p. 142) if =e is the amplification factor and the numerical value of actually obtained, an error of the form nmUnmU=nnnimmmUU eh; =1i with  a complex number and  real is considered. Then Equa-tion (2.15) leads to the following stability equation 210 32022410603 9052240027 807 142sin2sin=1 2816sin sin 4242sinsin sin4rr irrr rrrrirrrr2   (2.16) where 0 a typical value of , ; U1nmUnmU=0,1, ,m1N used for the linearization of the nonlinear terms, 00=U, 000=1  and =2h with 0,π2. Equation (2.16) is of the form =;,ABAB (2.17) with the set of the complex numbers, so the von Neumann necessary criterion for stability 1 will always be satisfied when .BA (2.18) Inequality (2.18) for =0 leads to 210 6010121 ,rr r   (2.19) which for =0 holds, while for =0 will be satis-fied when 1. (2.20) If =π2, inequality (2.18) leads to 251060390310116 2814 ,rr rrrrr    which subject to (2.20) holds. 2.3. The Modified Predictor-Corrector Scheme To avoid solving the nonlinear system (2.14) the follow-ing Modified Predictor-Corrector (MPC) scheme is pro-posed. 2.3.1. Predictor ˆtU is evaluated from the reccurence relation (2.12) replacing the matrix-exponential term with the following explicit second order rational approximant  22 21ˆ=02tIDDtas UU . (2.21) Then Equation (2.21) subject to Equation (2.11) using Copyright © 2011 SciRes. AJCM A. G. BRATSOS 155again the notations (2.7) - (2.8) leads to   2222ˆ = 2 .nnnnnn nnnn nttABtnABAABABABt   UUUUB (2.22) Let 21=ph, 2=2ph, 22 23=8ph, 224=2p, 245=2ph, 236=4ph, 27=ph4 and 238=4ph. Equation (2.22), when applied to the general mesh point of the grid G, gives  111 1211312211 1245211 262ˆ=2 464 2nn nnn nnmmmmmmmnnnnn nmmmmm mnnnnnnmmmmm mnnnnnmmmmmnnmmUUUpUUUpUUpUUUppUUU UUpU        11271111612 1 111 12811 1172 22 2 nnnmmmnnn nnmmm mmn nnnnnm mmmmmnnnnnmmmmmnn nnnnmmmm mmnnmmUUUpU UUpUU UUUpUU Up      1181 1 2; = 0,1,...,1.nnmmnnn nmmm mUUpU UUmNnmU (2.23) Stability Analysis Following again the Fourier method of analysing stabil-ity Equation (2.23) leads to the following stability equa-tion 220401 80242530 02670 6=14 2sin 16422sin sin42 sin24sin4ppppp ippp p ,  (2.24) which is of the form (2.17) with =1A. Then condition (2.18) for =0 leads to 222001112  . (2.25) which for =0 is obvious, while for =0 is satis-fied when the condition (2.20) holds. When =π2 condition (2.18) leads to 22 20002424841112hhh ,   (2.26) which again subject to (2.20) is always satisfied. 2.3.2. Corrector The corrector arises from Equation (2.13) as follows  222211 ˆ=21211 .212tDDtIDD  UUUt (2.27) Instead of the classical substitution of t  a m1nmUU the right-hand side of (2.27) by ˆtUodified pre-dictor-corrector method (MPC) was applied . The MPC method, which is explicit and is applied once, con-sists of considering (2.27) component-wise and using an updated component in the corrector vector as soon as it becomes available. Hence, in computing in, the cor-rected value 11nmU instead of the predicted value 11ˆnmU is used. The stability analysis of the corrector is given in Section 2.2.1. 3. Numerical Results For the linearization was given. Let the error at time level ; be 0=0,1,...,10=maxmNmUu=tn=1n,2,=0,1, ,1===maxeetLnnmu UmNm and ex the x-co- ordinate at which e occurs. Then e(2.2) - (2.3) denotes the error arising when using the boundary values (2.2) - (2.3), while analogous notations for the other boundary condi-tions are used. In all experiments the initial condition (1.2) was given by the value = ,0fxux=0.1h with u the theoretical solution (1.3). Experiments proved that the most accurate results are obtained for and 4=10 . For reasons of comparison with the correspond- ing works in [12,13,16,17] the same parameter values were used. 3.1. Problem  From the experiments the following are deduced: 1) when =0 (Table 1) using: i) the boundary values (2.2) - (2.3) the method intro-duced gives more accurate results for all time levels used than the corresponding results in  and marginally more accurate than those in [13,17], ii) the boundary condition (2.5) gives more accurate results than those in  and approximately equivalent to those in [13,17]. Copyright © 2011 SciRes. AJCM A. G. BRATSOS Copyright © 2011 SciRes. AJCM 156 From (i) - (ii) it is deduced that the boundary values (2.2) - (2.3) give more accurate results than the boundary condition (2.5). 2) when =0 (Table 2) using the boundary values (2.2) - (2.3) the method introduced has given -for =1 more accurate results for all time levels used than the corresponding in , and -for >1 results with marginally inferior accuracy to those in . In Figure 1(a) the solution u for 0,1t and 4410 ,10x is shown, while in Figure 1(b) the rele-vant solution U when 0,1x. Table 1. Problem . Comparisons of the proposed method for various values of x, t with α = 1, β = 1, γ = 0.001 and δ = 1 (h = 0.1, ℓ = 10–4). t x Exact e(2.2) - (2.3) e(2.5) e  e  e  0.1 0.5000187E−03 1.87406E−08 1.26463E−09 1.93715E−07 1.87406E−08 1.87405E−08 0.5 0.5000687E−03 1.87399E−08 1.97698E−08 1.93730E−07 1.87406E−08 1.87405E−08 0.05 0.9 0.5001187E−03 1.87250E−08 4.60177E−08 1.93745E−07 1.87406E−08 1.87405E−08 0.1 0.5000250E−03 3.74813E−08 6.39532E−09 3.87434E−07 3.74812E−08 3.74813E−08 0.5 0.5000750E−03 3.74736E−08 3.99558E−08 3.87464E−07 3.74812E−08 1.37481E−08 0.1 0.9 0.5001250E−03 3.74186E−08 7.66328E−08 3.87494E−07 3.74812E−08 3.74813E−08 0.1 0.5001374E−03 3.74814E−07 3.29223E−07 3.87501E−06 3.74812E−07 3.74812E−07 0.5 0.5001874E−03 3.72103E−07 3.79222E−07 3.87531E−06 3.74812E−07 3.74813E−07 1 0.9 0.5002374E−03 3.68427E−07 4.29222E−07 3.87561E−06 3.74812E−07 3.74813E−07 Table 2. Problem . Comparisons of the proposed method for various values of x, t and δ with α = 0, β = 1 and γ = 0.001 (h = 0.1, ℓ = 10–4). δ = 1 δ = 2 δ = 3 t x e(2.2) - (2.3) e  e(2.2) - (2.3) e  e(2.2) - (2.3) e  0.1 2.49875E−08 1.87465E−07 1.11763E−06 5.58901E−07 3.96731E−06 1.9841E−06 0.5 2.49875E−08 1.87486E−07 1.11750E−06 5.58836E−07 3.96652E−06 1.98371E−06 0.05 0.9 2.49874E−08 1.87508E−07 1.11737E−06 5.58772E−07 3.96572E−06 1.98331E−06 0.1 4.99750E−08 3.74934E−07 2.23526E−06 1.11779E−06 7.93462E−06 3.96811E−06 0.5 4.99750E−08 3.74977E−07 2.23500E−06 1.11766E−06 7.93304E−06 3.96731E−06 0.1 0.9 4.99749E−08 3.75019E−07 2.23474E−06 1.11753E−06 7.93144E−06 3.96652E−06 0.1 4.99750E−07 3.75002E−06 2.23526E−05 1.11754E−05 7.93462E−05 3.96632E−05 0.5 4.99749E−07 3.75044E−06 2.23500E−05 1.00741E−05 7.93303E−05 3.96553E−05 1 0.9 4.99749E−07 3.75086E−06 2.23474E−05 1.11728E−05 7.93143E−05 3.96473E−05 (a) (b) Figure 1. Problem  with δ = 1, α = 1, β = 1, γ = 0.001 when t  [0,1]: In (a) the surface shows u(x,t) for x  [–104, 104], while in (b) the numerical solution U whe n x  [0,1]. A. G. BRATSOS 157Table 3. Problem . Comparisons of the proposed method for various values of δ and γ when α = β = 1 (h = 0.1, ℓ = 10–4). t = 1 δ = 1 γ = 10–3 t = 0.5 δ = 2 γ = 10–2 t = 0.5 δ = 4 γ = 10–2 x e(2.2) - (2.3) e  x e(2.2) - (2.3) e  x e(2.2) - (2.3) e  0.1 3.74814E−07 3.74812E−07 0.1 3.89463E−05 2.75734E−04 0.1 5.69322E−05 1.08762E−030.5 3.72103E−07 3.74814E−07 0.3 3.89656E−05 2.75614E−04 0.3 5.69778E−05 1.08644E−030.9 3.68427E−07 3.74813E−07 0.5 3.89844E−05 2.75493E−04 0.5 5.70134E−05 1.08527E−03 Table 4. Problem . Boundary conditions (2.4 ) – (2.5). Comparisons of the propose d method for various values of t withα = 5 and δ = 1 (h = 0.1, ℓ = 10–4). γ = 10–3 γ = 10–4 γ = 10–5 t β Method e  Method e  e e  1 3.1570E−08 3.1616E−08 3.1584E−10 3.1630E−10 3.3410E−12 3.1632E−12 10 3.9684E−07 3.9742E−07 3.9702E−09 3.9760E−09 3.9704E−11 3.9762E−11 0.3 100 5.0291E−06 5.0365E−06 5.0316E−08 5.0389E−08 5.0318E−10 5.0392E−10 1 3.3393E−08 3.3394E−08 3.3408E−10 3.3409E−10 3.3410E−12 3.3411E−12 10 4.1976E−07 4.1977E−07 4.1995E−09 4.1996E−09 4.1997E−11 4.1998E−11 0.9 100 5.3165E−06 5.3166E−06 5.3221E−08 5.3223E−08 5.3224E−10 5.3225E−10 3.2. Problem  From Table 3 it is deduced that the method introduced using the boundary values (2.2) - (2.3) has given more accurate results for all time levels and parameters used than the relevant method in . 3.3. Problem  For reasons of comparison with the relevant results in  the boundary conditions (2.4) - (2.5) with 0=gt and 0,ut 1=1,gtut were used. From Table 4 it is deduced that the proposed method: -has given marginally more accurate results to those in  for all time levels and ,  used, -for fixed ,  and ● , the accuracy increases and U tends to identify with u at long time level as  is refined, ● , as  increases, the accuracy decreases. 4. Conclusions An implicit finite difference scheme based on fourth- order rational approximants to the matrix exponential term was proposed for the numerical solution of the Bur- gers-Huxley equation. The resulting nonlinear scheme was solved using an improved predictor-corrector method. The computational efficiency of the proposed method given in detail in Section 3 was tested by comparing the numerical results to selected ones in [12,13,16,17] using both appropriate boundary values and conditions. Con-clusions for the boundaries used were derived. Since the real world problems lead to the numerical solution of nonlinear equations or systems of equations, the introduced low cost and easy-to-handle method en-ables us to obtain accurate solutions. 5. References  A. L. Hodgkin and A. F. Huxley, “A Quantitative De-scription of Ion Currents and Its Applications to Conduc-tion and Excitation in Nerve Membranes,” Journal of Physiology, Vol. 117, No. 4, 1952, pp. 500-544.  J. Satsuma, “Topics in Soliton Theory and Exactly Solv-able Nonlinear Equations,” In: M. Ablowitz, B. Fuch-ssteiner and M. Kruskal, Eds., World Scientific, Singa-pore City, 1987.  X. Wang, “Nerve Propagation and Wall in Liquid Crys-tals,” Physics Letters, Vol. 112A, No. 8, 1985, pp. 402-406.  X. Y. Wang, Z. S. Zhu and Y. K. Lu, “Solitary Wave Solutions of the Generalised Burgers-Huxley Equation,” Journal of Physics A: Mathematical and General, Vol. 23, No. 3, 1990, pp. 271-274. doi:10.1088/0305-4470/23/3/011  A. G. Bratsos, “A Fourth-Order Numerical Scheme for Solving the Modified Burgers Equation,” Computers and Mathematics with Applications, Vol. 60, No. 5, 2010, pp. 1393-1400. doi:10.1016/j.camwa.2010.06.021  R. Fitzhugh, “Mathematical Models of Excitation and Propagation in Nerve,” In: H. P. Schwan. Ed., Biological Engineering, McGraw Hill, New York, 1969, pp. 1-85.  P. G. Estévez and P. R. Gordoa, “Painlevé Analysis of the generalized Burgers-Huxley Equation,” Journal of Phys-ics A: Mathematical and General, Vol. 23, No. 21, 1990, Copyright © 2011 SciRes. AJCM 158 A. G. BRATSOS pp. 4831-4837.  P. G. Estévez, “Non-Classical Symmetries and the Sin-gular Manifold Method: The Burgers and the Bur-gers-Huxley Equations,” Journal of Physics A: Mathe-matical and General, Vol. 27, No. 6, 1994, pp. 2113- 2127.  O. Yu. Yefimova and N. A. Kudryashov, “Exact Solu-tions of the Burgers-Huxley Equation,” Journal of Ap-plied Mathematics and Mechanics, Vol. 68, No. 3, 2004, pp. 413-420. doi:10.1016/S0021-8928(04)00055-3  X. Deng, “Travelling Wave Solutions for the Generalized Burgers-Huxley Equation,” Applied Mathematics and Computation, Vol. 204, No. 1, 2008, pp. 733-737. doi:10.1016/j.amc.2008.07.020  A. M. Wazwaz, “Analytic Study on Burgers, Fisher, Huxley Equations and Combined Forms of These Equa-tions,” Applied Mathematics and Computation, Vol. 195, No. 2, 2008, pp. 754-761. doi:10.1016/j.amc.2007.05.020  H. N. A. Ismail, K. Raslan and A. A. A. Rabboh, “Adomian Decomposition Method for Burger’s Huxley and Burger’s-Fisher Equations,” Applied Mathematics and Computation, Vol. 159, No. 1, 2004, pp. 291-301. doi:10.1016/j.amc.2003.10.050  I. Hashim, M. S. M. Noorani and M. R. S. Al-Hadidi, “Solving the Generalized Burgers-Huxley Equation Us-ing the Adomian Decomposition Method,” Mathematical and Computer Modelling, Vol. 43, No. 11-12, 2006, pp. 1404-1411. doi:10.1016/j.mcm.2005.08.017  M. Javidi, “A Numerical Solution of the Generalized Burger’s-Huxley Equation by Pseudospectral Method and Darvishi’s Preconditioning,” Applied Mathematics and Computation, Vol. 175, No. 1, 2006, pp. 1619-1628. doi:10.1016/j.amc.2005.09.009  M. Javidi, “A Numerical Solution of the Generalized Burgers-Huxley Equation by Spectral Collocation Method,” Applied Mathematics and Computation, Vol. 178, No. 2, 2006, pp. 338-344.  M. Javidi and A. Golbabai, “A New Domain Decomposi-tion Algorithm for Generalized Burgers-Huxley Equation Based on Chebyshev Polynomials and Preconditioning,” Chaos Soliton & Fractals, Vol. 39, 2009, pp. 849-857.  B. Batiha, M. S. M. Noorani and I. Hashim, “Application of Variational Iteration Method to the Generalized Bur-gers-Huxley Equation,” Chaos Solitons & Fractals, Vol. 36, No. 3, 2008, pp. 660-663. doi:10.1016/j.chaos.2006.06.080  A. Khattak, “A Computational Meshless Method for the Generalized Burgers-Huxley Equation,” Applied Mathe-matical Modelling, Vol. 33, No. 9, 2009, pp. 3718-3729. doi:10.1016/j.apm.2008.12.010  E. Babolian and J. Saeidian, “Analytic Approximate So-lutions to Burgers, Fisher, Huxley Equations and Two Combined Forms of These Equations,” Communications in Nonlinear Science and Numerical Simulation, Vol. 14, 2009, pp. 1984-1992. doi:10.1016/j.cnsns.2008.07.019  B. Fornberg, “A Practical Guide to Pseudospectral Meth-ods,” Cambridge University Press, Cambridge, 1998.  E. H. Twizell, “Computational Methods for Partial Dif-ferential Equations,” Ellis Horwood Limited, Chichester, 1984. Copyright © 2011 SciRes. AJCM