Research Article | Open Access
Development of Galerkin Method for Solving the Generalized Burger's-Huxley Equation
Numerical treatments for the generalized Burger's—Huxley GBH equation are presented. The treatments are based on cardinal Chebyshev and Legendre basis functions with Galerkin method. Gauss quadrature formula and El-gendi method are used to convert the problem into a system of ordinary differential equations. The numerical results are compared with the literatures to show efficiency of the proposed methods.
Hodgkin and Huxley  proposed a model in order to explain the ionic mechanisms underlying the initiation and propagation of action potentials in the squid giant axon. Satsuma investigated in 1986 that the generalized Burger’s-Huxley equation (GBH) shows a prototype model for describing the interaction between reaction mechanisms, convection effects, and diffusion transports . The most general form of GBH is  with the following initial and boundary conditions: where is a real parameter, , , and . If , then (1) becomes the generalized Burgers’ equation which is solved in . For , it changes to the generalized Huxley equation . Many researchers established numerical solutions of the GBH equation by using various numerical techniques. In  spectral collocation method is applied. Moreover, Darvishi’s preconditionings are employed to reduce round-off error in this method. In  a new domain decomposition algorithm based on Chebyshev polynomials and preconditioning is presented for solving GBH equation. In  up to tenth-order finite difference schemes are proposed to solve the GBH equation. In  Chebyshev spectral collocation with the domain decomposition is applied to find the numerical solution of the GBH equation.
Recently, in , a fourth-order finite difference scheme in a two-time level recurrence relation is proposed for the numerical solution of the GBH equation. The local discontinuous Galerkin method for Burger’s-Huxley equation has been studied in .
However, the various Galerkin algorithms have been applied in [12–14] for the numerical solutions of the ordinary differential equations. In this paper, the GBH equation is solved by nodal Galerkin methods. These methods have the advantages of both Galerkin and collocation methods. Collocation methods are derived from a strong form of the PDE and share many of the same advantages and disadvantages. Foremost is that they are easy to derive and to implement for a wide class of problems like constant coefficient, variable coefficient, and nonlinear problesms. Since collocation methods require the solution to satisfy the PDE at a set of grid points, they are naturally nodal approximations and will have aliasing errors even for constant coefficient problems. A main tradeoff is that there is little formal mathematical guidance on how to derive a stable approximation or how to implement boundary conditions. The latter makes it difficult to extend collocation methods to complex geometries or to systems of equations. Like finite difference methods, collocation methods are most easily applied to geometries that we can map onto a simple square or cube. Within those constraints, however, collocation methods will give spectrally accurate approximations.
Galerkin methods are derived from a weak form of the equations. They are less easily derived than collocation methods, but the formulation naturally leads to stable approximations and gives guidance on how to implement boundary conditions. Galerkin methods can be either nodal or modal. Modal approximations can be significantly more accurate than nodal approximations, depending on the problem. They are much harder to derive and more complex to implement, however, particularly for variable coefficient, nonlinear, or multidimensional problems.
Nodal Galerkin methods are intermediate between collocation and Galerkin methods. These methods start with the Galerkin formulation and the solutions are presented in a nodal form. Quadratures approximations are used to find the integrals that arise. The result is nodal methods which are significantly easier to implement than the Galerkin method and it can be extended to solve problems in complex geometries .
Furthermore, the presence of nonlinear term complicates the computation of the stiffness matrix [15, 16]. So nodal Galerkin methods are used and results are compared with a fourth-order finite difference scheme which implies to a nonlinear system  but our schemes imply to linear system which is easy to solve. Also, we compare our results with a new domain decomposition algorithm based on Chebyshev polynomials and preconditioning . In this study, the spectral collocation methods with the fourth-order Runge-Kutta method for time integration are used to solve the GBH equation. Moreover, preconditioning with the domain decomposition method is employed to reduce the round-off error in spectral collocation (pseudospectral) method. However, in our schemes we do not divide the domain into subdomains and we do not use any preconditioning to the resulted system. On the other hand, the domain decomposition method demanded small-step times to reach a good accuracy but in our methods we use big-step time and arrive to the same accuracy.
The remainder of this paper is organized as follows. In Section 2, we present Galerkin method with Chebyshev cardinal function as a basis function and we give the solution at the Chebyshev Gauss-Lobatto points. In Section 3, the Chebyshev Galerkin method with El-gendi quadrature is presented. In Section 4, El-gendi Legendre Galerkin method is described and we will use the Legendre cardinal function and the approximate solution will be presented at the Legendre Gauss-Lobatto points. In Section 5, numerical experiments are given and comparisons with the literatures to illustrate the efficiency of our methods are presented.
2. Gauss Chebyshev Galerkin (GCG) Method
In this section we explain the Gauss Chebyshev Galerkin method and illustrate how it is used to solve the problem (1) and (2) in case and . Let us now define some functional spaces. Let be a weight function on the interval and is the Banach space of the measurable functions such that The space is a Hilbert space for the following inner product: where Now, for any nonnegative integer the space is the space of all functions such that the derivatives of of order up to can be represented by functions in which is associated with the following norm: We also denote in particular the following space: Now, Let be any positive integer and the space of polynomials of degree at most ; we set We started by considering the approximation: where denotes the approximate value of , is the set of appropriate polynomials of degree , and is a set of coefficients. To have the approximation that satisfies the boundary conditions, we set . Now, the weighted Galerkin method  takes the form:
find Now we will take advantage of the flexibility given to us by the nodal representation of a polynomial. Since satisfies (10) for any function , it must satisfy the same condition for all linear combinations of the test functions: where are arbitrary coefficients. Since the test function is an th-order polynomial, we can write it in the following form: where is the cardinal Chebyshev polynomial and the nodal values are arbitrary, except that to ensure that satisfies the boundary conditions. We denote the approximation of the solution at the discrete grid points by , where (1) and (2) are enforced at the collocation points where Since the Chebyshev Gauss quadrature formula is given as follows: where ’s are given by In the Gauss Chebyshev Galerkin method, the trial function space coincides with the test function space which is a finite dimensional subspace of and is given as follows: whereas is given by for all , except and The grid points are called Chebyshev Gauss-Lobatto points, which are the extremal points of the Chebyshev polynomial . To get the nodal Galerkin approximation, we replace the integrals in (10) by Chebyshev Gauss-Lobatto quadrature.
Then the first discrete inner product becomes since , then the sum reduces to and the second inner product is If we rename the indices and , then we have where the first derivative of the cardinal functions at the points have the entries of the differentiation matrix : the second derivatives are In the same way, we can find the third and fourth term in By using (20), (22), (23), (24), and (25), we can write the discrete weak form in the following form: Since ’s are linearly independent, the coefficient of each must be zero, so Notice that the end points, and , are not included, since satisfying the boundary conditions. We specify the unknowns at those points by the boundary conditions. We complete the approximation (27) by using the approximation (9). So we have where The resulted system of ODEs has been solved by using fourth-order Runge-Kutta solver.
3. El-Gendi Chebyshev Galerkin (ECG) Method
In this method, the trial and test spaces are identical, so that we define for the space to be a vector space of functions such that all distributional derivatives of of order up to can be represented by functions in which is a Hilbert space for the inner product: Since the functions of are continuous up to the boundary by Sobolev imbedding theorem, it is meaningful. So, the solution subspace of is given as follows: which is a Hilbert space for the inner product defined as follows: where .
The weak forms of (1) and (2) are given by the following: find such that In this section, El-gendi formula has been used as follow: Let where are given by : In this case we use the following space: Now the discrete weak form is given as follows: find where Then the first term is given as follows: Similarly as above we can evaluated the rest of terms in the discrete weak form and we can write the resulting system after some manipulations as follows: where The resulted system of ODEs has been solved by using fourth-order Runge-Kutta solver.
4. El-Gendi Legendre Galerkin (ELG) Method
Analogous to the previous section we consider the Legendre cardinal function based on Legendre Gauss-Lobatto (LGL) nodes. El-gendi approximation will be used with a linear combination of the Legendre cardinal function as follows : where are the Legendre Gauss-Lobatto points and satisfies the condition: Now, the discrete weak form in the case of Legendre Galerkin method is as follows: find where Then the first discrete inner product becomes Since , then the sum reduces to and the second inner product is If we rename the indices and , then we have where is the first-order differentiation matrix that depends on Legendre polynomial at the LGL nodes and has the entries given by : where Also, has the following formula: where Then by using equations (45), (47), (49), (50), and (52) we can evaluate all the terms in the weak form () and we arrive to the following system: where As before, we use the fourth-order Runge-Kutta method to solve the resulted system (54).
Remark. In case of nonhomogenous boundary conditions the trial and test function space are the Sobolev spaces. So the boundary conditions are accounted naturally in the weak formulation. Moreover, the finite dimensional subspaces of the Sobolev spaces will be the whole space of polynomials. On the other hand, the methods which depend on numerical integration enforce the boundary condition explicitly or enforce the boundary points a particular linear combination of the approximate equation and the boundary condition; see [15, page 130-131] for more details.
5. Numerical Experiments
In this section we will give two examples and we will use MATLAB 7.0 software to obtain the numerical results. Consider the GBH equation: where , , , and are arbitrary constants and and are real numbers.
The initial condition is where The boundary conditions are where and the exact solution is The following error notations are defined: where and are the exact and approximate solutions, respectively. Since the Gauss-Lobatto points are in interval , the interval is mapped to by a linear mapping defined by where are the Gauss-Lobatto nodes.
Example 1. We consider the Burger’s-Huxley equation (1) with , , , , , , and . In Table 1, we compare the maximum error obtained by our method and domain decomposition algorithm in . As can be seen from this table, our methods are more accurate in time .
To show the solitary wave evolution with time, we expand the computation domain to and we plot the numerical solution and exact solution in Figure 1 for the values =1, , , and at . In Figure 2, we expand the computation domain to and plot the numerical solution and exact solution for the values , , and .
Example 2. In Table 2, the problem solved for various , and for , , , , , , and .
From Table 2 it is deduced that the proposed methods give more accurate results to those in  and the time step is not small as in . Also, it is noted that the accuracy decreased when increased and the accuracy increased when decreased.
Example 4. In Table 4, we will solve the same problem but with .