Abstract
This paper presents a new application of the homotopy analysis method (HAM) for solving evolution equations described in terms of nonlinear partial differential equations (PDEs). The new approach, termed bivariate spectral homotopy analysis method (BISHAM), is based on the use of bivariate Lagrange interpolation in the so-called rule of solution expression of the HAM algorithm. The applicability of the new approach has been demonstrated by application on several examples of nonlinear evolution PDEs, namely, Fisher’s, Burgers-Fisher’s, Burger-Huxley’s, and Fitzhugh-Nagumo’s equations. Comparison with known exact results from literature has been used to confirm accuracy and effectiveness of the proposed method.
1. Introduction
The study of nonlinear evolution partial differential equations (PDEs) is a vast area of research with well-developed and documented theories and applications in almost all areas of science and engineering. The PDEs are used to describe many complex nonlinear settings in applications such as vibration and wave propagation, fluid mechanics, plasma physics, quantum mechanics, nonlinear optics, solid state physics, chemical kinematics, physical chemistry, population dynamics, and many other areas of mathematical modelling. The development of both analytical and numerical methods for solving complicated highly nonlinear PDEs continues to be a fertile area of research geared towards enriching and deepening our understanding of these intriguing nonlinear problems.
The homotopy analysis method (HAM) has been widely discussed in the literature for solving both nonlinear ordinary and partial differential equations. A comprehensive exposition of the underlying concepts and applications of the HAM can be found in recently published books [1–4]. A unique feature of the HAM, which sets it apart from all perturbative and nonperturbative methods reported in the literature, is the flexibility to vary its embedded convergence controlling auxiliary parameters and functions. Previous studies have illustrated that a carefully selected choice of auxiliary linear operators result in significant improvement of the convergence and accuracy of the HAM [5–7]. Baxter et al. [5] examined multiple auxiliary linear operators to find the best operator that yields the best accuracy for the solution of the Cahn-Hilliard equation, a nonlinear partial differential equation. The linear operators in the study of Baxter et al. [5], as well as in numerous other HAM based studies for solving nonlinear PDEs, are conveniently chosen to guarantee that the HAM algorithm yields analytical series solutions. The limitation of these HAM based approaches that seek to obtain completely analytical results is that generating incremental terms of the HAM series solution becomes progressively cumbersome and the problem solving exercise becomes intractable eventually. This is particularly true when a nontrivial linear operator is used or required for optimal accuracy. Accordingly, an approach that can admit any form of linear operator, no matter how complex, is required in the HAM algorithm. However, complicated linear operators preclude the possibility of resolving the HAM series solutions analytically. Motsa et al. [8, 9] proposed a discrete version of the HAM that is based on Chebyshev spectral collocation approach for implementing the HAM algorithm which was otherwise impossible to solve analytically. This discrete variant of the HAM was called the spectral homotopy analysis method (SHAM) in [8, 9]. The SHAM has recently been extended to solve a nonlinear partial differential equation based problem of unsteady boundary layer flow caused by an impulsively stretching plate in [10]. Motsa [10] used a partial differential equation based auxiliary operator to improve on the ordinary derivative based linear operator approach used previously by Liao [11] to solve same problem. Motsa [10] concluded that when solving nonlinear PDEs, the use of PDE based linear operators leads to better results than the use of ODE based linear operators. In implementing the method, the decomposed equations were solved by applying the spectral method in the space variable and monomial series expansion in the time variable. This approach was found to work well for the unsteady boundary layer problem considered in [10] because the dimensionless time variable was defined in the range . However, series approaches of this kind are well known to have the capacity to resolve accurate solutions only in the region . Consequently, there is a need to develop variants of the SHAM that give solutions that are uniformly valid including regions where . More robust SHAM variations are needed for the solution of complex nonlinear PDEs that model important problems with wide applications in science, engineering, and other areas of applied mathematics.
The main objective of this work is to introduce a new variant of the spectral homotopy analysis method for solving nonlinear partial differential equations. The proposed method is developed by defining a rule of solution expression based on bivariate Lagrange interpolation. The homotopy analysis method algorithm is then applied to decompose the governing nonlinear PDEs into a sequence of linear PDEs. The resulting linear sequence of PDEs contains variable coefficients and is impossible to solve exactly. Consequently, the Chebyshev spectral collocation method is applied independently in the space and time independent variables. In view of the application of the combination of bivariate interpolation and spectral collocation differentiation, the new method is called bivariate interpolated spectral homotopy analysis method (BI-SHAM). The study presents a general BI-SHAM algorithm that can be used to solve second order nonlinear evolution equations. The applicability, accuracy, and reliability of the proposed BI-SHAM is confirmed by solving the Fisher, Burger-Fisher, Burger-Huxley, and Fitzhurg-Nagumo equations. The BI-SHAM results are compared against known exact solutions that have been reported in the scientific literature.
The remainder of the paper is organized as follows. In Section 2, we introduce the algorithm of the Bi-SHAM for a general nonlinear evolution PDE. Section 3 describes the application of the BI-SHAM in the problems that are selected for numerical experimentation. The numerical simulations and results are presented in Section 4. Finally, we conclude and describe future work in Section 5.
2. Bivariate Interpolated Spectral Homotopy Analysis Method (BI-SHAM)
In this section we introduce the Bivariate Interpolated Spectral Homotopy analysis Method (BI-SHAM) used to solve the governing nonlinear evolution PDEs. Without loss of generality, we consider nonlinear PDEs of the form where is the solution and is a nonlinear operator which contains all the spatial derivatives of . The solution procedure is based on the initial assumption that the solution can be approximated by a bivariate Lagrange interpolation polynomial of the form which interpolates independently at selected points in both the and directions defined as follows: The choice of grid points (3), which are called Chebyshev-Gauss-Lobatto points, ensures that there is a simple conversion of the continuous derivatives, in both space and time, to discrete derivatives at the grid points as will be discussed later. The functions are the characteristic Lagrange cardinal polynomials defined as follows: which obey the Kronecker delta equation. Consider The functions are defined in a similar manner.
To derive the HAM equations corresponding to the nonlinear equation (1), it is convenient to rewrite the governing equation in the form where the dot and primes denote the time and space derivatives, respectively, is a linear operator, and is a nonlinear operator. A crucial step in the implementation of the solution procedure is the evaluation of the time derivative at the grid points (). Before the derivative is applied, the given physical region, say , is converted to the region using the linear transformation . The values of the derivatives at the Chebyshev-Gauss-Lobatto points are computed as follows: where , are entries of the standard Chebyshev differentiation matrix of size (see, e.g., [12, 13]). By evaluating (6) at the grid points , we obtain
If the initial condition for (1) is given at (corresponding to ), we write (8) as follows: where is the known initial condition.
Equation (9) forms a system of coupled nonlinear ordinary differential equations with unknowns , . Below, we describe the spectral homotopy analysis method used to solve (9).
The algorithm of the HAM begins with the construction of the homotopy for a given linear operator defined as follows: where is an embedding parameter; denotes a nonzero convergence controlling auxiliary parameter; is the initial approximation of the solution of for . It should be emphasised that, in the specialized language of the HAM, the homotopy equation (10) is referred to as the zeroth order deformation equation. From (10) it can be noted that, as increases from 0 to 1, varies from the initial approximation to the solution of nonlinear equation (9). Expanding using Taylor series about gives where Thus, since and , we obtain The series (15) converges when the auxiliary parameter is carefully chosen. The functions appearing in series (15) are obtained as solutions of the so-called higher order deformations which are obtained by differentiating the zero-order deformation equation (10), times with respect to , then dividing by , and finally setting . This gives where The initial approximations are chosen in such a way that which, when using the definition of in (11), can be written as follows: Equation (20) to be solved for the initial approximations together with higher order deformation equations giving constitute sequence of linear ordinary differential equations and are solved using the Chebyshev spectral collocation method which is applied independently in the direction (where ) using Chebyshev-Gauss-Lobatto points. Consider defined as The derivatives with respect to are defined is terms of the Chebyshev differentiation matrix as where is the order of the derivative, () with being an Chebyshev derivative matrix, and the vector is defined as
Thus, substituting (23) in the equations that give the initial approximations (20), we obtain the following matrix system: where and is the identity matrix of size , the superscript denotes transpose, and the function is the coefficient of after the spectral method has been applied to the linear function . Solving (25) gives the initial approximation . To obtain the approximate solutions for (for ), the spectral collocation method, with discretisation in the direction, is applied in a similar manner to the higher order deformation equations (16). This gives the following matrix system: where and are as defined in (26) and In the above equation, is obtained by converting the continuous derivatives in to Chebyshev spectral derivatives.
3. Numerical Experiments
To demonstrate the applicability of the proposed Bi-SHAM algorithm as an appropriate tool for solving nonlinear partial differential equations, we apply the proposed algorithm to well-known nonlinear PDEs of the form (1) with exact solutions. In order to determine the level of accuracy of the BI-SHAM approximate solution, at a particular time level, in comparison with the exact solution we report maximum error which is defined by where is the solution obtained by (28) and is the exact solution at the time level .
Example 1. We consider Fisher’s equation as follows: subject to the initial condition and exact solution [14] where is a constant. Fisher’s equation represents a reactive-diffusive system and is encountered in chemical kinetics and population dynamics applications. For this example, the linear and nonlinear operator to be used in the Bi-SHAM algorithm are chosen as follows: Thus, using (17) we obtain The approximate solution at the is obtained by solving (28).
Example 2. We consider the generalized Burgers-Fisher’s equation [15] as follows: with initial condition and exact solution where , , and are parameters which, for illustration purposes, are chosen to be one in this paper. In this example, when , the linear and nonlinear operator to be used in the Bi-SHAM algorithm are chosen as follows: Thus, using (17) we obtain
Example 3. Consider the Fitzhurg-Nagumo equation as follows: with initial condition and exact solution [16] where is a parameter. In this example the linear and nonlinear operator to be used in the Bi-SHAM algorithm are chosen as follows: Thus, using (17) we obtain
Example 4. Consider the Burgers-Huxley’s equation as follows: where are constant parameters, is a positive integer (set to be in this study), and .
The exact solution subject to the initial condition, is reported in [17, 18] as where and . The general solution of the exact solution (48) was reported in [19, 20].
In this example the linear and nonlinear operator to be used in the Bi-SHAM algorithm are chosen as follows: Thus, using (17) we obtain
4. Results and Discussion
In this section we present the numerical solutions of the implementation of the BI-SHAM algorithm on the nonlinear evolution equations as described in the previous section. The number of collocation points in the space variable used to generate the results presented here was in all cases. Furthermore, unless otherwise specified, the order (total number of terms of the HAM series) of the HAM series was set to be . It was found that sufficient accuracy was achieved using these values in all numerical computations of the examples considered in this paper.
Using finite terms of the SHAM series we define th order approximation at the collocation points and (for and ) as follows:
Assuming that is the BI-SHAM approximate solution at the collocation (grid) points, the residual error is defined as where is defined as
The residual error is used in establishing the suitable convergence controlling parameter . A carefully selected is paramount in obtaining accurate and converging SHAM series solutions. The infinity norm of the residual error at a particular time level, defined as was used to identify the optimal value of that gives the best accuracy.
In Figures 1, 2, 3, and 4 we give illustrations of typical residual error curves that can be used to calculate the optimal value of for the Fisher, Burgers-Fisher, Fitzhugh-Nagumo, and Burgers-Huxley equations, respectively, when . The residual -curves are plotted using different orders of the BI-SHAM series. The optimal values of are chosen to be the clearly defined minimum of the residual curve. It can be seen from the figures that the optimal value lies in the range for Fisher’s equation, Fitzhugh-Nagumo’s equation, and Burgers-Fisher’s equation. In the case of Burgers-Huxley’s equation, it can be observed from Figure 4 that the optimal value is near −1.1. We also note that the residual error decreases with an increase in the order of the BI-SHAM series. This denotes convergence of the proposed method. It was also observed that convergence also improves with an increase in , the number of collocations used in the -variable. This observation is in accord with an earlier observation made in a related study by [21] where an interpolation based spectral homotopy analysis method was used to solve PDE based unsteady boundary layer flows.
In Tables 1, 2, 3, and 4 we give the maximum errors between the exact and BI-SHAM results (defined using (30)) for the Fisher, Burgers-Fisher, Fitzhugh-Nagumo, and Burgers-Huxley equations, respectively, at selected values of for different collocation points, , in the -variable. It is worth mentioning here that the results in all tables were computed on the space domain . To give a sense of the computational efficiency of the proposed method, the computational time taken to generate the results is also displayed in the tables. The results displayed in Tables 1, 2, 3, and 4 clearly show the accuracy of the proposed method. The accuracy is seen to improve with an increase in the number of collocation points . It is remarkable to note that accurate results with errors of order up to are obtained using very few collocation points in both the and variables , . This is a clear indication that the BI-SHAM is powerful method that is very appropriate in solving nonlinear PDEs of the type discussed in this investigation. We remark, also, that the BI-SHAM is computationally fast as the desired accurate results are generated in a fraction of a second in all the examples considered in this work.
5. Conclusion
This paper has presented a new variant of the spectral homotopy analysis method for solving general nonlinear evolution partial differential equations. The new method, called bivariate spectral homotopy analysis method (BISHAM), was developed from a combination of the homotopy analysis method algorithm with bivariate Lagrange interpolation and spectral collocation differentiation. The main goal of the current study was to assess the accuracy, applicability, and effectiveness of the proposed method in solving nonlinear partial differential equations. Numerical simulations were conducted on the Fisher equation, Burger-Fisher equation, Fitzhurg-Nagumo, and Burger-Huxley equations. This study has shown that the BISHAM gives very accurate results in a computationally efficient manner. Further evidence from this study is that the BISHAM gives solutions that are uniformly accurate and valid in large intervals of the governing space and time domains. The apparent success of the method can be attributed to the use of the nontrivial linear operators and the spectral collocation method for differentiating. This work contributes to the existing body of literature on homotopy analysis method based tools for solving complicated nonlinear partial differential equations. Further work needs to be done to establish whether the BISHAM can be equally successful in solving higher order nonlinear partial differential equations and coupled systems of two of more equations.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work is based on the research supported in part by the National Research Foundation of South Africa (Grant no. 85596).