Recent Advances in Solution Methods for Nonlinear Evolution Equations, Fluid Flow, and Heat and Mass Transfer
View this Special IssueResearch Article  Open Access
On the Optimal Auxiliary Linear Operator for the Spectral Homotopy Analysis Method Solution of Nonlinear Ordinary Differential Equations
Abstract
The purpose of this study is to identify the auxiliary linear operator that gives the best convergence and accuracy in the implementation of the spectral homotopy analysis method (SHAM) in the solution of nonlinear ordinary differential equations. The auxiliary linear operator is an essential element of the homotopy analysis method (HAM) algorithm that strongly influences the convergence of the method. In this work we introduce new procedures of defining the auxiliary linear operators and compare solutions generated using the new linear operators with solutions obtained using wellknown linear operators. The applicability and validity of the proposed linear operators is tested on four highly nonlinear ordinary differential equations with fluid mechanics applications that have recently been reported in the literature. The results from the study reveal that the new linear operators give better results than the previously used linear operators. The identification of the optimal linear operator will direct future research on further applications of HAMbased methods in solving complicated nonlinear differential equations.
1. Introduction
The homotopy analysis method (HAM) is an analytic method that has been widely used for solving highly nonlinear equations with applications in computational and applied mathematics, economics and finance, engineering, and many other areas of fundamental science. A systematic exposition of features of the HAM and its application can be found in recent books by Liao [1–3] and Vajravelu and Van Gorder [4]. A distinctive characteristic of the HAM that sets it apart from all other analytical methods is the presence of a convergencecontrolling parameter and the flexibility to select auxiliary functions and linear operators in order to guarantee convergence and improve accuracy of the approximate solutions. This makes the HAM suitable for solving many highly nonlinear problems including those that do not contain small or large embedded parameters.
A discrete numerical approach that uses the principles of the traditional HAM and combines them with the Chebyshev spectral collocation method was introduced by Motsa et al. [5, 6] and called spectral homotopy analysis method (SHAM). In the traditional HAM approach the options of initial approximations and auxiliary linear operators and functions are limited by the requirement that the obtained approximate solution must be a continuous analytical series solution. In the SHAM application, this restriction is removed and, consequently, the SHAM can accommodate infinite options of initial approximations, auxiliary linear operators, and convergencecontrolling auxiliary functions. This makes the SHAM more robust and capable of solving a wider range of complicated nonlinear equations than its analytical counterpart.
The effect of different auxiliary linear operators on the accuracy and convergence of the HAM has been studied by different authors in the past. Van Gorder and Vajravelu [7] presented different methods for selecting auxiliary linear operators and gave the necessary and sufficient conditions for the convergence of the HAM series solutions. The aim of the illustrative study of [7] was to demonstrate how different linear operators can be obtained. A comparative analysis of the accuracy and convergence of the approximate solutions obtained using their proposed linear operators was not considered. Shan and Chaolu [8] proposed a linear operator constructed by forming a differential equation using the guessed initial solution that satisfies the boundary conditions. The method proposed in [8] was demonstrated using the Thomas Fermi equation and a nonlinear heat transfer equation. Van Gorder [9] discussed methods for selecting auxiliary linear operators in solving the Painlevé equations, LaneEmden equation, and a third order equation that models a nonlinear stretching sheet in a fluid flow problem. It was noted that, for some auxiliary linear operators, the rate of convergence may be slow and for others the rate of convergence was rapid. All the studies reviewed so far, however, suffer from the fact that the linear operators have to be carefully selected in such a way that the solution method gives analytical solutions made up of simple functions such as polynomials, decaying exponential, sines, cosines, rational functions, and other elementary functions or products of such functions which are fairly easy to integrate. A particular method of selecting a linear operator such as the method of complete differential matching used in [7, 9] can only be used in very few nonlinear equations if the output solution is expected to be a simple series function. Despite being suggested in [7] as a choice of linear operator that, potentially, can give better convergence and accuracy, the method of complete differential matching has not been widely used in the HAM literature because it often leads to higher order decomposed equations that cannot be integrated analytically.
In seeking to identify auxiliary parameters and operators that give the best accuracy and convergence of the HAM algorithm, a more suitable study would seek to identify the optimal parameters from the infinitely many possibilities and not to search from the small set of options that guarantee exact solutions of the decomposed equations. The aim of this paper is to determine the auxiliary linear operator and initial approximation that give the optimal accuracy and convergence in the HAM solution of highly nonlinear boundary value problems defined in bounded domains. The study gives a systematic way of selecting initial approximations, auxiliary linear operators, and convergence controlling parameter. New methods of defining linear operators in the context of the discrete spectral method based version of the HAM (SHAM) are proposed. The accuracy and convergence of the SHAM algorithms derived using the proposed new linear operators are validated on the wellknown nonlinear differential equations that have been previously solved using the HAM/SHAM algorithms with the standard linear operators. In particular, we consider the DarcyBrinkmanForchheimer model [5, 10], JefferyHamel equation [6, 11–13], the laminar viscous flow in a semiporous channel subject to a transverse magnetic field [14–16], and the model of twodimensional viscous flow in a rectangular domain bounded by two moving porous walls [17–20].
The remainder of the paper is organized as follows. In Section 2, we give a brief description of the SHAM for general nonlinear ordinary differential equations. Section 3 introduces the auxiliary linear operators that are used in the SHAM algorithm. Section 4 describes the application of the SHAM with the proposed linear operators in the problems that are selected for numerical experimentation. The numerical simulations and results are presented in Section 5. Finally, we conclude and describe future work in Section 6.
2. Description of the Spectral Homotopy Analysis Method in a General Setting
In this section we give the basic description of the homotopy analysis method for the solution of a general nonlinear differential equation of order (where is a positive integer) given by where and are known functions of the independent variable and is an unknown function. The primes in (1) denote differentiation with respect to and . The function represents the nonlinear component of the governing equation. For illustrative purposes, we assume that (1) is to be solved in the domain subject to the separated boundary conditions where and are linear operators.
In the framework of the homotopy analysis method (HAM) [1–4], we define the following zeroth order deformation equations: where denotes an embedding parameter, is a kind of continuous mapping function of , and is the convergencecontrolling parameter. The nonlinear operator is defined from the governing equation (1) as with
By differentiating the zeroth order equations (3) times with respect to , setting , and finally dividing the resulting equations by , we obtain the following th order deformation equations: where
From the solutions of (6), the approximate solution for is determined as the series solution
A HAM solution is said to be of order if the above series is truncated at , that is, if
For nontrivial linear operators, the higher order equations (6) cannot be integrated using analytical means. Consequently, numerical approaches are employed. When the Chebyshev spectral collocation method is used to solve (6), the method is called spectral homotopy analysis method [3, 5, 6].
In using the SHAM, the initial guess is obtained simply as a solution of the linear part of the governing equation (1) subject to the underlying boundary conditions (2). That is, we solve
In essence, collocation approximates the solution using an interpolating polynomial of degree which satisfies the boundary conditions and the differential equations at all points, called the collocation points, , where . In the Chebyshev spectral collocation method, the collocation points are chosen to be the extrema of Chebyshev polynomials of degree on the interval defined as We use the transformation to map the interval to . The socalled differentiation matrix is used to approximate the derivatives of the unknown variables at the collocation points as the matrix vector product where and is the vector function at the collocation points. Higher order derivatives are obtained as powers of ; that is,
The matrix is of size and its entries are defined in [21, 22]. In the SHAM algorithm the continuous derivatives of the higher order deformation equations are replaced by the discrete Chebyshev differentiation matrices. As a result, the higher order deformation equations reduce to matrix equations and are solved using standard techniques for solving linear systems of equations.
3. Definition and Selection of Linear Auxiliary Operators
In this section we describe the different approaches used to define the auxiliary linear operators in the SHAM algorithm. In order to highlight the subtle differences between the variety of linear operators, it is convenient to express the nonlinear operator of the governing equation (1) using a sum formula of the derivative products as follows: where , are coefficients of the nonlinear terms containing as the highest derivative. In this way, (1) becomes
The selection of auxiliary linear operators in the application of the homotopy analysis method was described by VanGorder and Vajravelu [7] in a fairly general setting. In particular, three methods, namely, the method of highest order differential matching, linear partition matching, and complete differential matching, were defined in [7]. Below, we give a definition of the linear operator selection of [7] and present other options that can be implemented under the SHAM algorithm for any nonlinear differential equation that can be represented by (16).
(i) Method of the Highest Order Differential Matching. This approach defines the auxiliary linear operator using only the highest order derivative. This approach is the most widely used approach in the HAM solution of nonlinear differential equations defined in finite domains. With reference to (16), the linear operator is selected as (ii) Method of Linear Partition Matching. This approach defines the auxiliary linear operator to be the collection of all linear terms in the governing equation. With reference to (16), we set (iii) Method of Complete Differential Matching. In this approach, the collection of all linear operators and some nonlinear factors of the governing nonlinear equations are used to define the linear operator. The aim is to ensure that all the terms of the governing nonlinear differential equation contribute to the auxiliary linear operator selected. The following rules were defined in [7]:(a)in the case where we have a term which is the product of derivatives, we take the higher order derivative in the term;(b)if the term has a product of derivatives with functions of the unknown function, we again take the highest order derivative in the term;(c)in the case where we have a nonlinear expression in just the unknown function, we take the function itself.
Thus, applying the above rules to (16), we set (iv) Linear Partition Mapping after Transformation. This approach uses linear partition mapping after the transformation , where the function is carefully chosen to satisfy the underlying boundary conditions. In earlier versions of the SHAM [5, 6], it was suggested that be defined as a solution of (10) which is solved subject to the given boundary conditions. Substituting in the governing equation (1) to express the equation in terms of gives
Thus, the linear operator is set to be (v) Modified Complete Differential Matching Method. This is a new approach that is proposed as a modification of the method of complete differential matching of [7]. In the original approach [7], from the group of nonlinear terms of the governing equation, only the highest derivative is selected. In the proposed approach we want to ensure that all the terms of the nonlinear groups that make up the terms of the differential equation contribute to the linear operator. We define the following rules which are a modification of the rules set in [7].(a)In the case where we have a term which is the product of derivatives, we take the higher order derivative in the term and approximate each function in the remaining derivative product by . For example, if the nonlinear product is , we set , where is the solution of (10).(b)If the term has a product of derivatives with functions of the unknown function, we again take the highest order derivative in the term and approximate each function in the remaining derivative product by . For example, if the nonlinear product is , we set , where is the solution of (10).(c)In the case where we have a nonlinear expression in just the unknown function, we take the function itself and approximate the rest of the functions by . For example, if the nonlinear product is , we set , where is the solution of (10).
4. Numerical Experiments
In this section we illustrate how using a different linear operator can significantly improve convergence and accuracy of the SHAM. For illustration purposes we consider examples of nonlinear differential equations that have previously been solved in the published literature using the HAM/SHAM with standard linear operators. It is worth mentioning that in the previous HAMbased investigations the linear operators were chosen in such a way that(i)the obtained approximate solution is a continuous analytical expression;(ii)the HAM series solution conforms to a predefined rule of solution expression.The spectral method based approach of [5, 6] was aimed at removing the above restrictions by considering linear operators which are defined from part of the governing differential equation. In this case, it was found that the linearised deformation equations could not be solved exactly and hence the use of the spectral method. Below, we give a systematic procedure for defining the linear operators for various nondifferential equations selected for numerical experimentation.
4.1. Example 1: DarcyBrinkmanForchheimer Equation
We consider the following DarcyBrinkmanForchheimer equation that models the steady state pressure driven fully developed parallel flow through a horizontal channel that is filled with porous media: where is the dimensionless Forchheimer number and is the porous media shape parameter. This problem was previously solved using the SHAM with the auxiliary linear operator selected using the method of linear partition matching in [5]. The equivalent problem cast in cylindrical coordinates was subsequently solved in [10] where the linear operator was selected as the linear part of the governing equation excluding the linear terms with variable coefficients. For this reason, we remark that the linear operator used in [10] is not part of the class of linear operators considered in this study. Using the systematic choices described above, we set The linear operator is obtained by using the linear partition mapping after using the transformation . The resulting equation is The function is chosen as a function that satisfies the boundary conditions and the linear part of the (22). That is, must be a solution of Thus, the initial approximation is chosen to be The initial approximation used in the SHAM algorithms that use the linear operators , , and is obtained by solving the linear equations:
Thus, the initial approximations corresponding to the SHAM with , , and are given, respectively, by where .
The initial approximation to use in conjunction with is obtained as a solution of the differential equation formed from the linear part of the (24). That is, we solve This equation is a linear equation with variable coefficients. Thus, it is solved using the spectral method as described in the earlier section above. Similarly, the initial approximation to use with is obtained as a solution of the differential equation with given by (26).
With the nonlinear operator defined as the homogeneous part of the governing nonlinear equation in each case, we obtain the following higher order deformation equations: where
The nonlinear operator and higher order deformation equation corresponding to is defined as where and is the right hand side of (29).
4.2. Example 2: JefferyHamel Equation
Here, we consider the JefferyHamel equation that models the steady twodimensional flow of an incompressible conducting viscous fluid between two rigid plane walls that meet at an angle . The rigid walls are considered to be divergent if and convergent if . The governing equation is defined as subject to the boundary conditions where is the Reynolds number. This problem was investigated using the standard HAM with the highest order differential matching linear operator in [11, 12]. A related method, called optimal homotopy asymptotic method (OHAM), was used to solve the same problem in [13], again with . In Motsa et al. [6], the SHAM was used with the auxiliary linear operator defined using the transformed linear partition mapping method . In this example, we explore the auxiliary linear operators defined as
The function is chosen as a function that satisfies the boundary conditions and the linear part of (35). Thus, is
The linear operator is obtained by using the linear partition mapping after using the transformation . The resulting equation is subject to where The initial approximation used in the SHAM algorithms that employ the linear operators , , and is obtained by solving the linear equations: Thus, the initial approximations corresponding to the SHAM with , , and are given, respectively, by where . The initial approximation to use with is obtained as a solution of the differential equation formed from the linear part of (39). That is, we solve Since the above equation has variable coefficients, it is solved using the spectral method. Similarly, the initial approximation to use with is obtained as a solution of the differential equation with given by (38).
With the nonlinear operator defined as the homogeneous part of the governing nonlinear equation in each case, we obtain the following higher order deformation equations where
The nonlinear operator and higher order deformation equation corresponding to is defined as where
4.3. Example 3: Laminar Viscous Flow in a Semiporous Channel Subject to a Transverse Magnetic Field
In this section, we consider the problem of laminar viscous flow in a semiporous channel subject to a transverse magnetic field given by subject to the boundary conditions
This example was previously investigated using the standard HAM with the highest order differential matching linear operator in [14, 15]. In [16] a blend between the SHAM and a spectral method based quasilinearisation method was used to solve the same problem. In this investigation, we consider the following auxiliary linear operators: The initial approximation used in the SHAM algorithms that use the linear operators , , and is obtained by solving the linear equations: Thus, the initial approximations corresponding to the SHAM with and are given by We note that in this example is the second solution in (54). The initial approximations corresponding to and are obtained by solving (53) numerically using the spectral method. Similarly, is obtained by numerically solving subject to the boundary conditions where The higher order deformation equations are defined as where
The nonlinear operator and higher order deformation equation corresponding to is defined as where
4.4. Example 4: TwoDimensional Viscous Flow in a Rectangular Domain Bounded by Two Moving Porous Walls
In this section we consider the problem of twodimensional viscous flow in a rectangular domain bounded by two moving porous walls. The governing equations are given in [17, 23] as subject to the boundary conditions where is the nondimensional wall dilation rate defined to be positive for expansion and negative for contraction and is the permeation Reynolds number defined positive for injection and negative for suction through the walls. The nonlinear equation (62) was investigated using the standard HAM with the highest order differential matching linear operator in [17, 18]. Rashidi et al. [19] used the optimal homotopy asymptotic method (OHAM) to solve the same problem, again using as linear operator. The SHAM with transformed linear partition method was used in [20].
In this investigation, we consider the following auxiliary linear operators:
The initial approximations used in the SHAM algorithms that use the linear operators , , , and are obtained by solving the linear equations:
Thus, the initial approximation corresponding to is as follows: The initial approximations corresponding to , , and are obtained by numerically solving (65). Similarly, is obtained by numerically solving where
The higher order deformation equations are defined as where
The nonlinear operator and higher order deformation equation corresponding to is defined as where
5. Results and Discussion
The SHAM algorithm using the proposed auxiliary linear operators was applied to the several test problems presented in the last section. In this section we present the numerical results that highlight the difference in the accuracy of the approximate solutions and convergence of the different SHAM algorithms. All the simulations were conducted using collocation points. The accuracy and convergence of the SHAM approximate series solution depend on careful selection of the auxiliary parameter . The residual error was used to determine the optimal values of the auxiliary parameter. Using finite terms of the SHAM series we define th order approximation at the collocation points for as where . Given that is the SHAM approximate solution at the collocation points, the residual error at the collocation points is defined as where is the nonlinear operator defining the governing nonlinear equation (see (4)). We define the infinity norm of the residual error as The optimal that gives the best accuracy was identified to be the value of located at the minimum of the graph of the variation of the infinity norm against .
Figure 1 gives the typical residual error curves that are obtained using the 20th order SHAM approximate solution of the DarcyBrinkmanForchheimer equation (Example 1). The minimum of the curve gives the best possible residual error for each auxiliary linear operator used in the SHAM algorithm. A comparison of the minima of the residual error curves gives an indication of which auxiliary linear operator is likely to give better accuracy when used with the optimal auxiliary linear operator in the SHAM algorithm. It can be seen from Figure 1 that the best accuracy, at the 20th order SHAM approximation, is obtained when using the linear partition mapping after transformation (). The highest order differential mapping () gives the least accuracy. Figure 2 shows the variation of the norm of the residual error at the optimal value of as a function of the SHAM order. This graph gives an indication of how the SHAM series converges with an increase in the number of iterations (terms used in the series). Rapid convergence is determined by attaining the lowest possible residual error after few iterations. In this example it can be seen from Figure 2 that the highest order differential mapping linear operator is the slowest in terms of convergence, followed by the complete differential mapping . The operator that yields the fastest convergence in this example is which gives full convergence after only 10 iterations. It can also be seen from Figure 2 that , and all converge to the same level which is slightly lower than that of , but slightly higher than . This result suggests that the SHAM algorithm that applies gives the best accuracy and when using the accuracy is the lowest, in comparison with the other 4 choices of linear operators in this example.
Table 1 gives the minimum number of terms of the SHAM series that are required to attain full convergence of the SHAM algorithm using the optimal values of the convergence controlling parameter . The corresponding maximum residual errors are also tabulated. It can be seen from Table 1 that the least error is obtained when using and the other linear operators give more or less comparable residual error values. The SHAM with requires the fewest number terms of the SHAM series to attain maximum accuracy. The SHAM with comes second, followed by , , and in terms of the minimum SHAM order required to give best accuracy. This observation is in accord with the trends displayed in Figure 2. The table also indicates that the number of terms of the SHAM series required for the SHAM approximation to attain maximum accuracy increases with an increase in the value of the Forchheimer number parameter . It is worth noting from the governing equation (22) that the parameter multiplies the nonlinear component of the equation. Table 2 shows the effect of increasing the porous media shape parameter on the SHAM order required to attain the best accuracy. It can be seen from the table that the effect of on convergence and accuracy is similar to that of . In particular, it can be noted that as increases more SHAM series terms are required to attain the maximum accuracy in all linear operators.


Figure 3 shows the curves for the JefferyHamel problem (Example 2) corresponding to the different SHAM schemes generated using all the linear operators using 20 terms of the SHAM approximation series. It can be seen from Figure 3 that, after 20 iterations, would give the best accuracy and gives the least accuracy. This is in accord with the observation made earlier in Example 1. It can also be observed that the behaviour of the SHAM with is approximately the same as that of the SHAM with . This can be noted from the observation that the curve for lies on top of that for . In Figure 4 the variation of the residual error is shown as a function of the SHAM order. It can be observed that, in all test linear operators, the residual error decreases with the number of terms used in the SHAM series. This demonstrates that the SHAM approximation converges in all cases considered in this example. The convergence is fastest for followed by and is slowest when is used in the SHAM algorithm. It can also be observed that the results for , , , and converge to approximately the same level while converges to a level lower than that of the other linear operators. Thus, in addition to converging much faster, gives results that are slightly more accurate than those obtained using the other linear operators.
Table 3 presents the minimum SHAM order required to give maximum accuracy at the the optimal values of the convergence controlling parameter and the corresponding residual errors. The results are generated for increasing values of the Reynolds number (). It can be seen from Table 3 that the minimum order of the SHAM required gradually increases with an increase in for both and . There is no significant increase in the minimum SHAM order required in the case of , , and . The results from Table 3 also reveal that the required minimum SHAM order of and is comparable particularly for larger values of . This confirms the observation made in Figures 3 and 4. Thus, it can be inferred that the behaviour of the SHAM with higher order differential mapping () is comparable with that of SHAM with linear partition mapping () in the solution of the JefferyHamel equation.
