On the Optimal Auxiliary Linear Operator for the Spectral Homotopy Analysis Method Solution of Nonlinear Ordinary Differential Equations

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 well-known 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 HAM-based methods in solving complicated nonlinear differential equations.


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][2][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 convergence-controlling 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 convergence-controlling 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 2 Mathematical Problems in Engineering 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, Lane-Emden 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 well-known nonlinear differential equations that have been previously solved using the HAM/SHAM algorithms with the standard linear operators.In particular, we consider the Darcy-Brinkman-Forchheimer model [5,10], Jeffery-Hamel equation [6,[11][12][13], the laminar viscous flow in a semiporous channel subject to a transverse magnetic field [14][15][16], and the model of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls [17][18][19][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.

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  (0) ≡ ().
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][2][3][4], we define the following zeroth order deformation equations: where  ∈ [0, 1] denotes an embedding parameter, (; ) is a kind of continuous mapping function of (), and ℏ is the convergence-controlling parameter.The nonlinear operator N is defined from the governing equation (1) as with By differentiating the zeroth order equations (3)  times with respect to , setting  = 0, 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  = 0, 1, . . ., .In the Chebyshev spectral collocation method, the collocation points are chosen to be the extrema of Chebyshev polynomials of degree   on the interval −1 ≤  ≤ 1 defined as We use the transformation  = ( − )( + 1)/ is the vector function at the collocation points.Higher order derivatives are obtained as powers of D; that is, The matrix  is of size (+1)×(+1) 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.

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  , (), ( = 0, 1, . . ., ) 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 Van-Gorder 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 () = ()+ 0 (), where the function  0 () is carefully chosen to satisfy the underlying boundary conditions.In earlier versions of the SHAM [5,6], it was suggested that  0 () 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  0 ().For example, if the nonlinear product is   ()  ()  (), we set   0 ()  0 ()  (), where  0 () 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  0 ().For example, if the nonlinear product is  2 ()  (), we set  2 0 ()  (), where  0 () 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  0 ().

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 HAM-based 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.

Example 1: Darcy-Brinkman-Forchheimer Equation.
We consider the following Darcy-Brinkman-Forchheimer 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 L 2 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 L 4 is obtained by using the linear partition mapping after using the transformation () = () + ỹ0 ().The resulting equation is The function ỹ0 () is chosen as a function that satisfies the boundary conditions and the linear part of the (22).That is, ỹ0 () 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 L 1 , L 2 , and L 3 is obtained by solving the linear equations: Thus, the initial approximations corresponding to the SHAM with L 1 , L 2 , and L 3 are given, respectively, by where The initial approximation  0 () to use in conjunction with L 4 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 L 5 is obtained as a solution of the differential equation with ỹ0 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 Mathematical Problems in Engineering The nonlinear operator and higher order deformation equation corresponding to  = 4 is defined as where and  1 () is the right hand side of (29).

Example 2: Jeffery-Hamel Equation.
Here, we consider the Jeffery-Hamel equation that models the steady twodimensional flow of an incompressible conducting viscous fluid between two rigid plane walls that meet at an angle 2.
The rigid walls are considered to be divergent if  > 0 and convergent if  < 0. The governing equation is defined as subject to the boundary conditions where Re is the Reynolds number.This problem was investigated using the standard HAM with the highest order differential matching linear operator L 1 () in [11,12].A related method, called optimal homotopy asymptotic method (OHAM), was used to solve the same problem in [13], again with L 1 ().In Motsa et al. [6], the SHAM was used with the auxiliary linear operator defined using the transformed linear partition mapping method L 4 ().In this example, we explore the auxiliary linear operators defined as The function ỹ0 () is chosen as a function that satisfies the boundary conditions and the linear part of (35).Thus, ỹ0 () is The linear operator L 4 is obtained by using the linear partition mapping after using the transformation () = () + ỹ0 ().The resulting equation is subject to where The initial approximation used in the SHAM algorithms that employ the linear operators L 1 , L 2 , and L 3 is obtained by solving the linear equations: Thus, the initial approximations corresponding to the SHAM with L 1 , L 2 , and L 3 are given, respectively, by where  1 = √ 4 2 + 2 Re.The initial approximation  0 () to use with L 4 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 L 5 is obtained as a solution of the differential equation with ỹ0 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  = 4 is defined as where (49)

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 L 1 () 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 L 1 , L 2 , and L 3 is obtained by solving the linear equations: L  ( 0 ) = 0,  0 (0) = 0,   0 (0) = 0,  0 (1) = 1,   0 (1) = 0, for  = 1, 2, 3.
(53) Thus, the initial approximations corresponding to the SHAM with L 1 and L 2 are given by We note that in this example ỹ0 is the second solution in (54).The initial approximations corresponding to L 3 and L 5 are obtained by solving (53) numerically using the spectral method.Similarly,  0 () 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  = 4 is defined as where (61)

Example 4: Two-Dimensional Viscous Flow in a Rectangular Domain Bounded by Two Moving Porous Walls.
In this section we consider the problem of two-dimensional 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 Re 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 L 1 () in [17,18].Rashidi et al. [19] used the optimal homotopy asymptotic method (OHAM) to solve the same problem, again using L 1 () as linear operator.The SHAM with transformed linear partition method L 4 () was used in [20].

Mathematical Problems in Engineering
In this investigation, we consider the following auxiliary linear operators: The initial approximations used in the SHAM algorithms that use the linear operators L 1 , L 2 , L 3 , and L 5 are obtained by solving the linear equations: (65) Thus, the initial approximation  0 () corresponding to L 1 is as follows: The initial approximations corresponding to L 2 , L 3 , and L 5 are obtained by numerically solving (65).Similarly,  0 () is obtained by numerically solving where The higher order deformation equations are defined as where (70) The nonlinear operator and higher order deformation equation corresponding to  = 4 is defined as where (72)

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  = 20 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  = 0, 1, 2, . . .,  as where Given that Y is the SHAM approximate solution at the collocation points, the residual error at the collocation points is defined as where N 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 Darcy-Brinkman-Forchheimer 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 (L 4 ).The highest order differential mapping (L 1 ) 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 L 1 is the slowest in terms of convergence, followed by the complete differential mapping L 3 .The operator that yields the fastest convergence in this example is L 4 which gives full convergence after only 10 iterations.It can also be seen from Figure 2 that L 2 , L 3 , and L 5 all converge to the same level which is slightly lower than that of L 1 , but slightly higher than L 4 .This result suggests that the SHAM algorithm that applies L 4 gives the best accuracy and when using L 1 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 L 1 and the other linear operators give more or less comparable residual error values.The SHAM with L 4 requires the fewest number terms of the SHAM series to attain maximum accuracy.The SHAM with L 5 comes second, followed by L 2 , L 3 , and L 1 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 Jeffery-Hamel 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, L 4 would give the best accuracy and L 1 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 L 1 is approximately the same as that of the SHAM with L 2 .This can be noted from the observation that the ℏ-curve for L 1 lies on top of that for L 2 .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 L 4 followed by L 5 and is slowest when L 1 is used in the SHAM algorithm.It can also be observed that the results for L 1 , L 2 , L 3 , and L 5 converge to approximately the same level while L 4 converges to a level lower than that of the other linear operators.Thus, in addition to converging much faster, L 4 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 (Re).It can be seen from Table 3 that the minimum order of the SHAM required gradually increases with an increase in Re for both L 1 and L 2 .There is no significant increase in the minimum SHAM order required in the case of L 3 , L 4 , and L 5 .The results from  Table 3 also reveal that the required minimum SHAM order of L 1 and L 2 is comparable particularly for larger values of Re.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 (L 1 ) is comparable with that of SHAM with linear partition mapping (L 2 ) in the solution of the Jeffery-Hamel equation.
Figure 5 shows the ℏ-curves corresponding to the example of laminar viscous flow in a semiporous channel subject to a transverse magnetic field (Example 3).Again, it can be seen from Figure 5 that the residual error is lowest when the SHAM is used with L 4 and largest when used with L 1 which, as was observed in the case of Jeffery-Hamel flow, is comparable with L 2 .The same observation can be made in Figure 6 where the residual error is plotted against the SHAM order.From Figure 6, it can also be seen that L 4 gives the best convergence followed by L 5 and the slowest convergence is obtained when using L 1 which is comparable to L 2 .Table 4 shows effect of increasing the Reynolds number (Re) on the minimum SHAM orders and the corresponding residual errors.It can be observed, from Table 4, that the minimum SHAM orders increase with an increase in Re when the SHAM is used with L 1 , L 2 , and L 3 .It can also be seen that there is no significant increase in the minimum SHAM orders when using L 4 and L 5 .Furthermore, the residual error obtained using L 4 is lower than that of the other linear operators.This is in line with similar observations made earlier.
Table 5 shows the effect of increasing the Hartmann number (Ha) on the minimum SHAM orders and the corresponding residual errors in Example 3. In this case, it can be observed that the minimum SHAM orders increase with an increase in Ha only in the case of L 1 .In the case of the other linear operators the minimum SHAM orders do not increase with an increase in Ha.This means that convergence of the SHAM series is independent of the value of Ha when the SHAM is used with L 2 , L 3 , L 4 , and L 5 .It can also be noted from Table 5 that the residual error is smallest when L 4 is used.If can be inferred that the SHAM is most accurate when used with L 4 .
Figure 7 shows the ℏ-curves corresponding to the example of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls (Example 4).It can be noted from Figure 7 that the residual error is the lowest when the SHAM is used with L 4 and the largest when used with L 1 after 20 iterations.Figure 8 illustrates the variation of the residual error against the SHAM order.It can be seen from Figure 8 that L 4 gives the best convergence followed by L 5 and the slowest convergence is obtained when using L 1 .The linear operators L 1 , L 2 , L 3 , and L 5 converge to the same level of residual error which is higher than that of L 4 .Again, it can be inferred that L 4 gives superior accuracy compared to the other linear operators.
Table 6 shows the effect of increasing the Reynolds number (Re) on the minimum SHAM orders and corresponding residual errors in Example 4. It can be observed that the minimum SHAM orders increase with an increase in Re when the SHAM is used with L 1 , L 2 , and L 3 .In the case of L 4 and L 5 it seems possible that an increase in Re has no strong influence in the convergence speed of the SHAM.It can also be noted from Table 6 that the residual error is the smallest when L 4 is used.This is in accord with the trends observed in all the other examples considered in this work.Table 7 shows the effect of the nondimensional wall dilation rate  on the minimum SHAM order and corresponding residual errors.It can be observed from Table 7 that an increase in  results in a corresponding significant increase in the minimum SHAM order for L 1 .The minimum SHAM orders are not significantly influenced by  when the SHAM is used with L 2 , L 3 , L 4 , and L 5 .Furthermore, the best accuracy is obtained when using L 4 .

Conclusion
This study gives a systematic way of choosing initial approximations and linear operators that can be employed in the spectral homotopy analysis method (SHAM) algorithm.
The main aim of the study was to determine the linear operator definition that gives approximate solutions with optimal convergence and accuracy.A comparison between the accuracy and convergence of SHAM algorithms using commonly used linear operators and newly proposed linear operators was conducted.Simulations were conducted on selected nonlinear boundary value problems defined on bounded domains.One of the more significant findings to emerge from this study is that the linear operator defined by the highest order differential matching, which is the most commonly used linear operator for solving boundary value problems defined on bounded domains, results in the slowest converging and least accurate SHAM algorithm.
The SHAM algorithm that uses a linear operator defined by the method of linear partition matching was found to perform slightly better than the SHAM with the highest order differential matching in some cases.However, in other cases the methods were found to be comparable.It was also shown that the SHAM with linear operator defined by the method of complete differential matching performed much better than the SHAM with linear partition mapping or the highest order differential matching.The study proposed a modification of the complete differential method approach which, unlike the original approach, does not ignore the nonlinear factors of the derivative terms in the definition of the linear operator.This new definition of linear operator was found to perform better than the original complete differential method approach.The most optimal linear operator was found to be one that uses the linear partition mapping after a transformation that expresses the unknown function as a sum of a new unknown plus some initial guess that satisfies the underlying boundary conditions.The current findings add substantially to our understanding of how homotopy analysis based methods can be modified to yield approximate solutions with optimum accuracy and convergence.Further work needs to be done to investigate the application of the linear operators proposed in this study in solving nonlinear partial differential equations.It would also be very interesting to extend the current study to systems of nonlinear differential equations.

1 Figure 6 :
Figure 6: Variation of residual error with SHAM order for Example 3 when Re = 10 and Ha = 1.