The nonlinear Harry-Dym equation is  ,
with the initial approximation,
where a and b are constants.
Equation (1) is one of the most important soliton equations i.e., a soliton is deemed as solitary, travelling wave pulse solution of nonlinear partial differential equations (NPDEs)  e.g., the Korteweg-de Vries (KdV) equation (see Refs.   ). Moreover, according to Drazin and Johnson  , a soliton can interact strongly with other solitons and retain its identity. Among the efforts and achievements due to the behavior of the Harry-Dym equation, Popowicz  generalized this equation by the Lax and Hamiltonian formulations. He also found that generalization of the Harry-Dym Lax operator is only limited to two operators. As reported by Marciniak and Błaszak  , if a classical Stäckel system is applied to construct the coupled Harry-Dym (cHD) hierarchy, both nonlocal and purely differential parts of hierarchies are obtained. Using the Riemann- Liouville derivative, Huang and Zhdanov  recently investigated the time fractional nonlinear Harry-Dym equation with the Lie symmetry and obtained some invariant solutions explicitly. Brunelli and da Costa  studied a large class of nonlocal equations and charges for the Harry-Dym hierarchy. Tian and Liu  presented a new non-trivial supersymmetric Harry-Dym equation using the supersymmetric reciprocal transformation. Ma  constructed an extended Harry-Dym hierarchy using a sequence of Lax equations with the new time variables. However, this extended Harry-Dym equation which is based on the adjoint eigenfunctions, have not been taken into account by Wu et al.  . In other words, Wu et al.  proposed a 2+1 dimensional Harry-Dym hierarchy and found that this extended scheme is Lax integrable. It should be noted that more details on the Harry-Dym equation are available in Refs.   .
This paper is organized as follows. Section 2 gives an overview of the HAM. In Section 3, the HAM is applied to solve the nonlinear Harry-Dym equation with the initial approximation given in Equation (2). Section 4 is devoted to the results and discussion. Concluding remarks are given in Section 5.
2. Homotopy Analysis Method (HAM)
The essence of the HAM  -  is based on the approximation of PDEs along with a convenient way to guarantee convergence of the series solution. This issue can be outlined by the following example.
Consider the following nonlinear differential equation,
in which A is a nonlinear operator, x and t are independent variables and is an unknown function. On account of a continuous mapping from , the homotopy embedding parameter p varies from 0 to 1; it means that varies from the initial approximation to the exact solution . This mapping is constructed by the zeroth-order deformation equation which is defined as,
in which L is a linear operator, is an auxiliary parameter and is an auxiliary function. Expanding in the Taylor series (see Appendix A) with respect to p gives,
If , , and are properly chosen, the power series (6) converges at , then one has,
which must be one of the solutions of Equation (3), as proved by Liao  . Now, it is of particular interest to show that Equation (7) provides a relationship between the initial approximation and the exact solution so that:
1) Solution of Equation (4) is valid for .
2) The deformation derivative exists for all values of m.
Differentiating the zeroth-order deformation Equation (4) m times with respect to p, then setting and dividing them by , gives the so-called mth-order deformation equation as,
It should be noted that for is governed by the linear Equation (8) with the linear boundary conditions that come from the original problem, which can be easily solved using some software programs such as MAPLE or MATHEMATICA.
Then can HAM be expeditiously used in NPDEs? There have been a number of NPDEs which can be solved through the HAM. To address this subject, Nassar et al.  developed the HAM to derive an approximate solution of the Poisson-Boltzmann equation for semiconductor devices. The solutions of the fractional Swift Hohenberg equation can be found in the work of Vishal et al.  . They could report influence of the real bifurcation parameter on probability density function and then concluded that the proper values of auxiliary and homotopy parameters are needed. However, their model was not capable of taking into account the two-dimensional domain for fractional Brownian motion. Hetmaniok et al.  took advantage of the HAM to solve linear and nonlinear integral equations of the second kind. They also used the auxiliary parameter to show influence of this parameter on the convergence rate. In a book, current advances of the HAM to solve strongly nonlinear problems were edited by Liao  . Ren et al.  introduced an approximate analytical solution to the Gross- Pitaevskii equation while a nonlinear Schrödinger equation was utilized to simulate Bose-Einstein condensates trapped in a harmonic potential. Martin  proposed an analytical algorithm to solve an integral-differential equation of the transport theory for stationary case using HAM. The solution of nonlinear boundary value problems (NBVPs) based on the Chebyshev operational matrices was investigated by Shaban et al.  . They applied the Tau method to convert a set of algebraic equations so that the solution can be obtained iteratively. Motsa et al.  utilized the spectral-homotopy analysis method (SHAM) to solve the Darcy-Brinkman-Forchheimer equation and found the admissible range of auxiliary parameter for this problem. Hesameddini and Latifizadeh  obtained solutions for the first and second types of the Painlevé equations with 8, 10 and 12 iterations. Wang  developed an approximate analytical solution of relativistic Toda lattice system and concluded that using HAM may be very effective for solving differential difference equations (DDEs). In other words, some papers are devoted to compare the HAM with the other analytical methods e.g., optimal homotopy asymptotic method (OHAM), differential transformation method (DTM), homotopy perturbation method (HPM), general series expansion method, harmonic balance method (HBM) which can be found in Refs.      .
3. HAM Formulation for the Nonlinear Harry-Dym Equation
As mentioned above, the Equation (1) can be solved through the HAM. To this end, the linear operator takes the form of,
with the property,
in which d is an integral constant. Now the nonlinear operator can be expressed in terms of Equation (1) as,
Using the above definition, with the assumption , the zeroth-order deformation can be given as,
For and , it is respectively written as,
It should be noted that the discretized form of above equation is illustrated in Appendix B. So, solution of the
Now, according to Equation (18), the values of for can be obtained as,
To continue this procedure for , can be achieved and the series solution is completely computed.
4. Results and Discussion
4.1. Comparison and Validation
To validate the present iterative approach, results are compared with the exact solution of the Harry-Dym equation provided by Mokhtari  which can be expressed as,
where c is a constant. Mokhtari  also showed that the variational iteration method (VIM) is not able to obtain implicit solution of the Harry-Dym equation. However, one should point out that Hereman et al.  proposed implicit solution of this equation via a direct method and also performed transformations between the Harry-Dym and KdV equations explicitly. In particular, according to the fractional Harry-Dym equation,
where is a parameter describing order of the fractional derivative, Figure 1 depicts variations of versus x for the Harry-Dym equation calculated by different methods. As shown in this figure, variation of predicted by the HAM is in good agreement with those provided by Mokhtari  . However, there is a relative error between the results which can be calculated by,
where and are the exact and approximate solutions of Equation (1), respectively. Due to and  , it is noted that the maximum error does not exceed 1.098%, as shown in Figure 2. The root-mean- square error (RMSE) between the exact solution and the present iterative
Figure 1. Variation of versus x for 7th-order approximation of the nonlinear Harry-Dym equation with the HAM and compared with those reported in Refs.      .
Figure 2. Depiction of discrepancy (i.e., error) between the present iterative method and exact solution due to Equation (22).
approach, i.e.,  , is also 0.509%. Apart from the
exact solution, other observations may be made from Figure 1; e.g., it is assumed that for , the series solution will be converged at the results presented in   , and in general, the HAM consists of the HPM and Adomian decomposition method (ADM). In fact, the two exhibit similar pattern for different values of the auxiliary parameter. Therefore, it is particularly important to account for the interaction between the auxiliary parameter and initial approximation before using the HAM even for the cases which NPDEs have been taken into account. This type of scheme has been presented by Turkyilmazoglu  as an illustrative comparison between the HAM and HPM. In addition, using 7th- order approximation of the HAM instead of ADM presented in Ref.  is suggested largely due to its less time-consuming and also is able to obtain valid values of the auxiliary parameter. However, it should be noted here that although we utilize in our solution, the convergence can be observed when the objective is to replace Equation (21) by Equation (1) for .
Another comparison is also performed here to illustrate the accuracy of present iterative approach when the effect of nonlinearity has been taken into account. In connection with this matter, Rawashdeh  utilized the fractional reduced differential transformation method (FRDTM) to solve the nonlinear Harry-Dym equation while the maximum error between the HAM and FRDTM for is 11.723%. In addition, He found that predicted by the FRDTM is on average 1.18 less than that predicted by exact solution and is greatly affected by the number of iterations performed. Based on this comparison, it is seen that Equation (21) without would also fit into HAM, namely, Caputo sense  , which is not correspondent to define the fractional order initial conditions. Herein, the Caputo’s formula is defined as  ,
It should be noted that Kumar and Singh  applied the Laplace transform on the both sides of Equation (21) and solved the governing equation through a homotopy perturbation transform method (HPTM), but their solution could not capture especially for the cases with large values of auxiliary parameter. There exists a discrepancy between the exact solution and those obtained through HPTM because of decomposing the nonlinear term in He’s polynomial (i.e., Equation (17)). To clarify, the He’s polynomial has been also used in some previous studies   . It is to be noted that Guo and Mei  studied time- fractional Boussinesq-type equation by the VIM and Khan and Wu  investigated the advection problems through HPTM.
4.2. Universal Graphs for Convergence Region of the Presented Iterative Approach
The fact that Equation (18) is not convergent for all values of the auxiliary parameter, lifts some particular restrictions on choosing a proper value of (e.g., by trial and error  ). Above all, making a good selection of h is important not only with the convergence, but also with the basis functions   . It should be noted that other perturbation methods require the need for small parameters, which may commonly diverge series solutions with slight changes (see Ref.  ). Indeed, at large values of the auxiliary parameter, convergence may lead to only few terms of series solution. It is shown that based on the zeroth- order deformation Equation (4), one may expect the flexibility and freedom to select L, despite the fact that the order of linear operator will be different from the original nonlinear problem (see Ref.  ). As said in Ref. 
“∙∙∙ we are free to enhance the convergence region and convergence rate of a series solution via an appropriate choice of the auxiliary parameter even for fixed choice of the initial approximation, auxiliary linear operator and auxiliary function.”
―R.A. Van Gorder and K. Vajravelu 
Due to Figure 1 that varies almost linearly versus x, Figure 3 shows the variation of for some different auxiliary parameters. It is seen that indicates the ascending behavior while suggests the descending behavior of the present iterative approach. Although the slopes of these lines are slightly different from each other, the lines become parallel by increasing or decreasing the auxiliary parameter in the neighbourhood of . It should be noted that this fact can result in almost linear variation of . In theory, at the mth-order of approximation, the square residual error can be defined as  ,
It should be noted that as decreases to zero more rapidly, the convergence rate for corresponding series solution would become faster  . In this regard, Yabushita et al.  introduced an analytical approximation to projectile motion with the quadratic resistance law by the HAM. They derived governing equations using a power series which was accurate enough for large angle of projection. Therefore, with a minimum value of , the optimal values of unknown auxiliary parameters were obtained.
Figure 3. Universal graphs to show the effect of the auxiliary parameter on the convergence of present iterative approach.
It should be noted that by increasing t, will be decreased slightly, as it is shown in Figure 4(a) and Figure 4(c). It is mainly due to the fact that this parameter plays an important role in higher-order governing equations (e.g., Equation (19c)) (see Refs.     ). It is to be noted that a similar conclusion for increasing t can also be drawn through a similar procedure. Odibat  reached a pragmatic approach for the convergence of HAM, and to date, no counterexample exits in the open literature for this conclusion, neither for NPDEs nor for PDEs. Hence, in such graphs we increase t from 1 to its a neighborhood which can exhibit this deviation for a better presentation. However, it is to be noted that this parameter is small in most of the analytical approximations (e.g., )    .
According to the higher-order equations (e.g., Equation (19c)), if t will be increased, the maximum value of will be increased for (see Figure 4(b) and Figure 4(d)). It is noted that almost drops to zero when the auxiliary parameter reaches to zero as well. It is seen that neglecting the effect of t may lead to underestimation of versus h. However, if this parameter exceeds its , which can be named as a minimum value of t, will be fully negative.
It is noteworthy to mention that although many techniques (e.g., HPM, ADM, FRDTM, HPTM) can be useful for predicting analytical approximations of the nonlinear Harry-Dym Equation (1), the more accurate results can be obtained
Figure 4. Universal graphs for predicting the threshold of .
through the HAM. On the other hand, since the HAM is employed to approximate NPDEs analytically, it is important to note that it enables us to construct a mapping of initial approximation to the exact solution of nonlinear governing equation.
5. Concluding Remarks
In this paper, the nonlinear Harry-Dym equation, which describes a couple of nonlinearity and dispersion, has been solved iteratively. It has been shown that Equation (1), which can be transformed to the KdV equation, is an integrable NPDE. Furthermore, it has been shown that the HAM is an efficient method which provides a convenient way to guarantee convergence of the series solutions. The main inferences that can be drawn from this work are summarized as follows:
・ The HAM will enable us to make an informed choice for auxiliary linear operator and initial approximation which can be employed to obtain the series solution.
・ Equation (18) indicates one important role of the auxiliary parameter which was introduced to construct the zeroth-order deformation Equation (4).
・ Convergence of the power series (6) largely depends on the range of auxiliary parameter.
・ In case of , the HPM is a special case of HAM, as proved by Liao  .
・ By using , influence of the auxiliary function, with the exception of some particular recursive problems (see Ref.  ) will be neglected.
・ Presence of affects considerably when the fractional Harry- Dym equation is utilized.
We start with the case z, h and ( ) as complex numbers and all singular points of , respectively. If is analytic at , one has,
By setting in Equation (A1), we obtain the classical Taylor series (see Ref.  ).
The discretized form of Equation (17) for solving Equation (1) can be expressed as,