Research Article  Open Access
Development of Galerkin Method for Solving the Generalized Burger'sHuxley Equation
Abstract
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 Elgendi 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.
1. Introduction
Hodgkin and Huxley [1] 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â€™sHuxley equation (GBH) shows a prototype model for describing the interaction between reaction mechanisms, convection effects, and diffusion transports [2]. The most general form of GBH is [3] 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 [4]. For , it changes to the generalized Huxley equation [5]. Many researchers established numerical solutions of the GBH equation by using various numerical techniques. In [6] spectral collocation method is applied. Moreover, Darvishiâ€™s preconditionings are employed to reduce roundoff error in this method. In [7] a new domain decomposition algorithm based on Chebyshev polynomials and preconditioning is presented for solving GBH equation. In [8] up to tenthorder finite difference schemes are proposed to solve the GBH equation. In [9] Chebyshev spectral collocation with the domain decomposition is applied to find the numerical solution of the GBH equation.
Recently, in [10], a fourthorder finite difference scheme in a twotime level recurrence relation is proposed for the numerical solution of the GBH equation. The local discontinuous Galerkin method for Burgerâ€™sHuxley equation has been studied in [11].
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 [15].
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 fourthorder finite difference scheme which implies to a nonlinear system [5] 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 [7]. In this study, the spectral collocation methods with the fourthorder RungeKutta method for time integration are used to solve the GBH equation. Moreover, preconditioning with the domain decomposition method is employed to reduce the roundoff 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 smallstep times to reach a good accuracy but in our methods we use bigstep 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 GaussLobatto points. In Section 3, the Chebyshev Galerkin method with Elgendi quadrature is presented. In Section 4, Elgendi Legendre Galerkin method is described and we will use the Legendre cardinal function and the approximate solution will be presented at the Legendre GaussLobatto 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 [17] 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 thorder 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 GaussLobatto points, which are the extremal points of the Chebyshev polynomial . To get the nodal Galerkin approximation, we replace the integrals in (10) by Chebyshev GaussLobatto 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 [18]: 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 fourthorder RungeKutta solver.
3. ElGendi 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, Elgendi formula has been used as follow: Let where are given by [19]: 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 fourthorder RungeKutta solver.
4. ElGendi Legendre Galerkin (ELG) Method
Analogous to the previous section we consider the Legendre cardinal function based on Legendre GaussLobatto (LGL) nodes. Elgendi approximation will be used with a linear combination of the Legendre cardinal function as follows [20]: where are the Legendre GaussLobatto 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 firstorder differentiation matrix that depends on Legendre polynomial at the LGL nodes and has the entries given by [21]: 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 fourthorder RungeKutta 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 130131] 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 GaussLobatto points are in interval , the interval is mapped to by a linear mapping defined by where are the GaussLobatto nodes.
Example 1. We consider the Burgerâ€™sHuxley equation (1) with , , , , , , and . In Table 1, we compare the maximum error obtained by our method and domain decomposition algorithm in [7]. 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 .
(a)
(b)
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 [7] and the time step is not small as in [7]. Also, it is noted that the accuracy decreased when increased and the accuracy increased when decreased.
Example 3. In Table 3, the problem solved for various and for , , , , , , , and . In Table 3, we show maximum error and compare with the corresponding results in [10].

Example 4. In Table 4, we will solve the same problem but with .
