Explicit Solutions to Large Deformation of Cantilever Beams by Improved Homotopy Analysis Method I: Rotation Angle

: An improved homotopy analysis method (IHAM) is proposed to solve the nonlinear differential equation, especially for the case when nonlinearity is strong in this paper. As an application, the method was used to derive explicit solutions to the rotation angle of a cantilever beam under point load at the free end. Compared with the traditional homotopy method, the derivation includes two steps. A new nonlinear differential equation is ﬁrstly constructed based on the current nonlinear differential equation of the rotation angle and the auxiliary quadratic nonlinear differential equation. In the second step, a high-order non-linear iterative homotopy differential equation is established based on the new non-linear differential equation and the auxiliary linear differential equation. The analytical solution to the rotation angle is then derived by solving this high-order homotopy equation. In addition, the convergence range can be extended signiﬁcantly by the homotopy–P á de approximation. Compared with the traditional homotopy analysis method, the current improved method not only speeds up the convergence of the solution, but also enlarges the convergence range. For the large deﬂection problem of beams, the new solution for the rotation angle is more approachable to the engineering designers than the implicit exact solution by the Euler–Bernoulli law. It should have signiﬁcant practical applications in the design of long bridges or high-rise buildings to minimize the potential error due to the extreme length of the beam-like structures.


Introduction
In the past two decades, large-deformation beam structures have been widely applied in long-span bridges [1,2], high-rise buildings [3,4], aerospace, and multi-flexible robots. With the rapid development of very high-performance steel and sophisticated bridge structures, the increasing number of bridge constructions across rivers and seas are developing trendily towards longer, higher, and lighter structures. The span or height of the structure is increasing continually, and the structural configuration is becoming more and more sophisticated to cater for more stringently esthetic and functional demands. On the other hand, available analytical solutions to this special form of structure do not predict its performance explicitly with its implicit integral, which raises potential uncertainties in the deformations, especially when extreme large span is concerned. Therefore, it is imperative to establish more sophisticated theories and methods for the high structural nonlinearity in beam-like structures with extreme lengths.
The perturbation method [5][6][7][8] is one of the most prevailing analytic tools for nonlinear problems. It has played significant roles in the development of modern science and engineering as it footprints itself in revealing many important properties and interesting phenomena of non-linear problems. The perturbation technique is based on the existence of a small parameter, which is called a perturbation quantity. Though the perturbation quantity is a cornerstone of perturbation techniques, it also poses restrictions on the practical application of the perturbation techniques. It is not possible that each non-linear problem contains such a perturbation quantity, wherein the existence of perturbation quantities may not necessarily ensue a satisfactory result. For example, both the straightforward perturbation method and the singular perturbation technique fail to formulate effective theoretical drag solutions to a sphere in a uniform stream [9,10]. This is mainly because that, like other analytic methods, perturbation techniques cannot effectively control the convergence of approximation series and adjust convergence regions when necessary.
The unfavorable dependence of perturbation techniques on small parameters might be circumvented by introducing an artificial small parameter. In 1892, in consideration of the equation where A(t) is a time periodic matrix, Lyapunov [11] introduced an artificial parameter m. Thus, Equation (1) becomes The power series expansion is then calculated over a range of the value of m for the solutions. In many cases, Lyapunov proved that the series converge when m = 1, and thus, the final expression can be identified with m = 1. The above approach is called Lyapunov's artificial small parameter method [11]. The further development of this approach was registered by the δ-expansion method by Karmishin et al. [12]. On the other hand, both the artificial parameter method and the δ-expansion method fail to provide a principal to identify the domain and specify the artificial parameter or δ. Similar to perturbation techniques, both the artificial small parameter and the δ-expansion methods do not control effectively the convergence of approximation series and adjust convergence regions when necessary.
To circumvent the aforementioned improficiencies, the idea of an artificial parameter is generalized by homotopy [13]. The homotopy is a fundamental concept in topology [14] that can be traced back to Jules Henri Poincaré (1854~1912), a French mathematician. Based on homotopy, two methods have been developed. One of these derivants, the homotopy continuation method, dates back to the 1930s [15][16][17][18] and is a global convergent numerical method mainly for nonlinear algebraic equations. The other is the homotopy analysis method (HAM) proposed in the 1990s by Shijun Liao [19,20], which is an analytic approximation method with guarantees convergence mainly for nonlinear differential equations. The homotopy analysis method has been successfully applied to various types of ordinary differential equations and partial differential equations.
One of the successful applications of the homotopy analysis method [21,22] was fulfilled by introducing an auxiliary parameter ћ to construct a new general form of homotopy. Unlike other analytic techniques, this homotopy analysis method always serves a family of analytic results at any given order of approximation. Generally, the homotopy analysis method has the following advantages: It is valid even if a given non-linear problem does not contain any small perturbation parameters at all; ii.
It controls effectively the convergence of approximation series and adjusts convergence regions when necessary; iii.
It efficiently approximates a non-linear problem by choosing different sets of base functions.
The large deflection problem of beams is governed by the Euler-Bernoulli law. To derive elastic large deflections, the equivalent load, deflection, and longitudinal displacement can be expressed by analytic elliptic integrals with end point rotation angles as parameters [40]. For the cantilever beam under concentrated force, the solution to this elastic problem contains elliptic integrals, which is not explicit or straightforward for practical application [41][42][43]. The homotopy analysis method (HAM) was used to obtain the solution for the large deformation of a cantilever beam made of axially functionally graded material [44] and a nonlinear beam subjected to a coplanar terminal load consisting of a moment, an axial compressive force, and a transverse force [45]. Closed-form series solutions for the static deflection of anisotropic composite beams resting on elastic foundations were obtained by both the homotopy analysis method (HAM) and iterative HAM (iHAM) [46]. The iterative homotopy analysis method (iHAM) was used to obtain analytical solutions for the arbitrary large deflection of geometrically exact beams subjected to distributed and tip loads based on follower and conservative loading scenarios [47]. Wang et al. [48] derived an explicit solution to the large deformation of a cantilever beam under point load at the free end with HAM. However, their homotopy analysis method had the problem of slow convergence because it used iterative solutions of linear differential equations to approximate the exact solution to strong non-linear differential equations. In addition, explicit solutions to vertical and horizontal displacements were not given.
In this paper, an improved homotopy analysis method (IHAM) is proposed to solve strong nonlinear differential equation, and as an application, it is used to obtain an explicit solution to the rotation angle for the large deformation of a cantilever beam under point load at the free end. In another parallel paper, this improved homotopy analysis method is applied to derive the explicit solutions for vertical and horizontal displacements of large deformations of cantilever beams.

Formulations
The proposed homotopy analysis method (IHAM) is essentially different from the HAM in view of construction methods.
HAM constructs a high-order nonlinear iterative homotopy differential equation using selected linear differential equation and the original non-linear differential equation, whereas IHAM formulates a high-order nonlinear iterative homotopy differential equation with the selected linear differential equation, the selected simple non-linear differential equation, and the original non-linear differential equation. Table 1 lists the formulation steps of IHAM and HAM for comparison. Table 1. Improved homotopy analysis method vs. homotopy analysis method.
To approximate the exact solution, differences also rise from their respective approaching paths. HAM approximates the exact solution of the original non-linear differential equation by iterating with linear differential equation. Whereas IHAM uses linear differential equation and simple non-linear differential equation to iterate in approximating the exact solution to the original non-linear differential equation. Thus, the convergence rate of IHAM is much higher than that of HAM. IHAM not only accelerates the convergence rate, but also enlarges the convergence range in solving strong nonlinear problems.
In this study, a boundary value problem is considered for second order nonlinear differential equations, as can be seen in the following form Here f (w) is a derivable function with the initial values of f w (0) = c 0 , f w (0) = c 1 , and 1 2! f w (0) = c 2 . In compliance with the boundary conditions, we have It is worth noting that the prerequisite of any small or large perturbation parameters in Equation (3) is dismissed. Thus, the proposed approach is general.

Zero-Order Deformation Equation
For comparison, nonlinear boundary value Equations (3) and (4) are solved firstly by HAM. When w(ξ) is expanded into a power series of elementary functions in terms of ξ, we have where a k is elementary coefficient. Equation (8) serves as the approximation solution to the current problem. With reference to the boundary condition in Equations (4) and (8), leads to which will be specified as an initial guess solution for w(ξ). The homotopy analysis method is based on assumptive continuous mapping of w(ξ) → Ψ(ξ, q) . When the embedded variable q increases from 0 to 1, Ψ(ξ, q) changes from the initial guess solution w 0 (ξ) to the exact solution w(ξ). According to the Equations (3) and (8), the below auxiliary linear operator is selected It satisfies L(C 1 + C 2 ξ) = 0 where C 1 and C 2 are coefficients. From Equation (3), the following nonlinear operator is specified Denoting a non-zero auxiliary parameter as ћ, a more general form of homotopy is then constructed as Setting [Ψ(ξ, q)] = 0, a family of equations are derived as follows Equation (14) is subject to the boundary conditions where q ∈ [0, 1] is an embedded parameter. When q = 0, Equations (14) and (15) yield When q = 1, Equations (14) and (15) are exactly the same as Equations (3) and (4) Therefore, as the embedded parameter q increases from 0 to 1, Ψ(ξ, q) evolves from the initial guess w 0 (ξ) to the solution w(ξ). For brevity, Equations (14) and (15) are defined as zero-order deformation equations.
Assume that Ψ(ξ, q) is analytic in the domain of q ∈ [0, 1], so that the deformation derivatives exist: Ψ(ξ, q) can then be expanded in the Maclaurin series of q as follows: where It is worth mentioning that the auxiliary parameter ћ in Equation (14) determines the convergence regions. Assuming that ћ is properly chosen, so that all of these Maclaurin series are convergent at q = 1. Thus, at q = 1 Equation (17) gives Equation (21) bridges the initial guess solution w 0 (ξ) to the exact solution w(ξ). The result for the nth-order approximation is given by:

High-Order Deformation Equation
Define a vector w n = {w 0 (ξ), w 1 (ξ), w 2 (ξ), . . . , w n (ξ)} Differentiating Equations (15) and (16) n times with respect to q and then setting q = 0 and dividing them by n! gives the governing equation of w n (ξ), in conjunction with the boundary conditions and The right-hand side of Equation (24) can be derived by symbolic calculation software such as Maple, MATLAB, etc. Thereafter, the linear high-order deformation Equations (24) and (25) are solved.

Construction of a New Nonlinear Homotopy Differential Equation
Using the original differential Equation (3) and the selected differential Equation (7), the nonlinear operator is defined as: The improved homotopy analysis method replaces the nonlinear operator N in Equation (12) with the new nonlinear operator N in Equation (28).
When ε = 0, the new nonlinear operator N in Equation (28) degenerates to the nonlinear operator N in Equation (12).
When q = 1, the new nonlinear differential equations N = 0 becomes the original differential Equation (3) for N = 0.
When q = 0, the nonlinear differential equations N = 0 is the simple differential Equation (7) for N = 0.

Construction of High-Order Homotopy Equation
With the new nonlinear operator N in Equation (28) and the selected linear operator L in Equation (10), the homotopy operator is defined as follows, The improved homotopy analysis method replaces the homotopy operator H in Equation (13) with the new homotopy operator H in Equation (29).
When ε = 0, the new homotopy operator H in Equation (29) degenerates to the homotopy operator H in Equation (13).
With definition of a vector a new homotopy equation of high order is similarly constructed in conjunction with the boundary conditions and The improved homotopy analysis method replaces the homotopy Equation (14) with the new homotopy Equation (31) for high-order problem.
To solve Equations (31) and (32), the result of nth-order approximation are given by

Convergence Theorem
Theorem 1 (Convergence Theorem). If the series converges, where w n (ξ) satisfies the high-order deformation Equations (31) and (32), and the definitions in Equations (33) and (34) are true, it is the solution to Equations (3) and (4).
The above assumption is then proved as below. Letting being the convergent series, a necessary condition for the series to converge is: From Equations (31), (34) and (38), we have Generally, Ψ(ξ, q) does not satisfy the original nonlinear Equation (3). Letting represents the residual error of Equation (3), then δ(ξ, q) = 0 corresponds to the exact solution of Equation (3). According to the above definition, the McLaughlin series of residual error δ(ξ; q) with respect to the embedded variable q is: According to Equation (40), when q = 1, Equation (42) becomes: According to the definition of δ(ξ; q), when q = 1, Equation (43) is the exact solution to the original Equation (3). Therefore, as long as the below series: converges, there exists a solution to the original Equation (3). (36) converges, where w n (ξ) satisfies the high-order deformation Equations (31) and (32), and the definitions in Equations (33) and (34) are true, there is:

Theorem 2. If the series in Equation
This theorem is self-explanatory with reference to Equation (40). According to Theorems 1 and 2, only the initial guessing solution w 0 (ξ), auxiliary parameter ћ, auxiliary linear operator L, and auxiliary nonlinear operator N to ensure the convergence of series in Equation (36) should be derived.

Problem Description
We consider the large deformation of a cantilever beam under point load at the free end as shown in Figure 1. The bending equation of a uniform cross-section beam with a large deformation is: [40,41] dθ It complies to the boundary conditions: where s is the arc-coordinate of the neutral axis of the beam, x is the horizontal coordinate from the fixed end, L stands for the length of the beam, F denotes the point load at the free end, EI specifies the bending stiffness of the beam, θ represents the rotation angle of cross-section of the beam, and l specifies the unknown horizontal distance of two ends. Differentiating the equation with respect to S and then using the dimensionless variables ξ = S/L, the original Equation (45) becomes in compliance to boundary conditions of The prime in Equations (47) and (48) denotes differentiation with respect to ξ, and The rotation angle of cross-section plane at the free end is denoted by θ b = θ(1). For infinitesimal deformation, the below linear equation will be sufficient which is subjected to the boundary conditions The corresponding solution is: which gives the linear result of ration angle at the free end for ξ = 1

Zero Order Deformation Equation
The nonlinear boundary value Equations (47) and (48) are then solved by IHAM to derive a explicit formula for θ b . To start with, θ(ξ) is expanded into the power series of ξ: Equation (52) is then identified as an initial guess. Accordingly, the initial guess of θ(ξ) is expressed as: For simplicity, the below simplest auxiliary nonlinear differential equation is chosen which is subjected to the boundary conditions where ε is the correction coefficient of the convergent region, and ε ∈ [0, 1]. When ε = 0, Equation (56) becomes the linear Equation (50). When ε = 1, Equation (56) becomes: In fact, the approximating second-order Taylor expansion for the cosine function in Equation (47) can be written as: According to the original Equation (47) and the equivalent Equation (56), we define one nonlinear operator as: where q ∈ [0, 1] is an embedded parameter, ψ(ξ, q) is a function dependent on ξ and q, ε is a control parameter for the convergent region, and ε ∈ [0, 1]. When q = 0, there is: Being an auxiliary linear operator, and ћ ∈ [−1, 0) denotes a nonzero auxiliary parameter (convergence-control parameter), a homotopy is then constructed as: When q = 0, we have If q = 1, Equation (61) becomes: Thus, by enforcing We have a family of equations: which are subjected to the boundary conditions of: The prime in Equation (63) denotes differentiation with respect to ξ. When q = 0, Equation (62) incurs: When q = 1, Equations (62) and (63) are equivalent to the original Equations (47) and (48), provided that: Thus, as q increases from 0 to 1, ψ(ξ, q) evolves from the known initial guess θ 0 (ξ) to the solution θ(ξ) in Equations (47) and (48).
By expanding ψ(ξ; q) into a Taylor series of the embedded parameter q and using Equation (64), we have: where Assuming that ћ and ε are properly chosen so that the series Equation (66) are convergent at q = 1, according to Equation (65), we have the series: And the governing equations of θ n (ξ) can then be deduced from the zero-order deformation Equations (62) and (63).

High-Order Deformation Equations
Substituting Equation (66) into Equation (62) and differentiating Equation (62) n times with respect to the embedded parameter q, dividing it by n!, and then setting q = 0, we have the nth-order deformation equations as follows L[θ n (ξ) − χ n θ n−1 (ξ)] = ћR n (θ 0 , θ 1 , · · · , θ n−1 ) Equation (68) is subjected to the boundary conditions: where and Combining Equations (59) and (70), we have: The right-hand side of R n can then be calculated with symbolic calculation software such as Maple, MATLAB and so on. Thereafter, the linear high-order deformation Equations (68) and (69) are solved.
The solutions by IHAM contain the auxiliary parameter ε about nonlinear operator N 0 and the auxiliary parameter ћ about linear operator L, which control the convergence region and the rate of the IHAM solution series.

Calculation Results
For simplicity, the rotation angle with dimension one is defined as: It is important to ensure a series to be convergent in a sufficiently large region. Generally, the convergence region and the rate of the series are governed by the basic function approaching the solution. Unlike traditional analytic methods, the improved homotopy analysis method accommodates variable basic functions to approximate the solution of a given nonlinear problem. Therefore, the series of solutions that converge in the whole region of independent variables of physical significance can be obtained with IHAM. It is worth noting that the form of solution expression is crucial to the convergence of the problem since it specifies the initial guess solution θ 0 (ξ) and auxiliary linear operator L. It would be helpful to put forward the so-called solution expression principle: the initial guess solution θ 0 (ξ) and the solution θ 1 (ξ), θ 2 (ξ), . . . of the high-order deformation equation should not violate the solution expression.
Provided that the initial guess solution θ 0 (ξ) and the auxiliary linear operator L are chosen, the homotopy analysis method still affords flexibility in identifying the value of the auxiliary parameter ћ. Superior to all the traditional analytical methods, HAM always offers a series of solutions with auxiliary parameter ћ. The auxiliary parameter ћ dominates the convergence region and rate of the series of solutions. By selecting an appropriate ћ value, the convergence region and the rate of the series of solutions can be effectively controlled.
Similar to the auxiliary parameter ћ, even if the initial guessing solution θ 0 (ξ), auxiliary linear operator L and auxiliary parameter ћ are specified, again IHAM is flexible with the value of auxiliary parameter ε. Different from the original HAM, the IHAM gets a series of solutions with the variable auxiliary parameter ε, which also affects the convergence region and rate of the series of solutions. By selecting the appropriate ε value, the convergence region and the rate of the series of solutions will be effectively adjusted.
Therefore, superior to the traditional analytical methods, the improved homotopy analysis method efficiently controls the convergence region and the rate of the series of solutions by adjusting the auxiliary parameters ћ and ε.

Exact Solution to Rotation Angle
From Equation (47), the transcendental equation solved θ b is as follows where is the first type of elliptic integral and K(µ) is the first complete elliptic integral. Therein According to Equations (53) and (73), the linear solution to the rotation angle at the free end with dimension 1 is: It can be verified that the convergent region of the linear solution (78) is as follows if the relative error ∆ ≤ 1% is required.

Effect of Auxiliary Parameter ε on Convergence
It is expected that the improved homotopy analysis method derives a family of series solutions with the variable auxiliary parameter ε. The efficiency of the IHAM lies in the appropriate identification for the value of ε to ensure a convergent solution promptly in a sufficiently large region.
The parameters in nonlinear problems could contain physical connotations, such as the angular frequency of nonlinear vibration, the wall friction of viscous flow, etc. These parameters are usually correlated to the auxiliary parameter ε. Therefore, the curve of these physical quantities can be drawn with respect to ε by regarding the auxiliary parameter ε as an independent variable. For example, below Θ b is a physical quantity If Θ b is a function of, ε, Θ b − ε can be identified. According to Theorem 1 or Theorem 2, all convergent series of Θ b given by different values of ε converge to the exact solution. If the solutions are unique, they should all converge to the same value. Therefore, there is a horizontal segment in the Θ b − ε curve to imply the converge region of ε represented by R ε . For simplicity, such a curve is called a ε curve, and the corresponding region R ε is the effective region of ε. Obviously, if the value of ε falls in the valid region, the series solution will converge without fail. Even though Θ b does not make any physical sense, the corresponding ε curve can still be attempted. The more ε curves that are drawn, the more specific value of ε will be identified. If the initial guess solution θ 0 (ξ), auxiliary linear operator L and the auxiliary parameter ћ are given, the effective region of ε for different physical quantities is usually almost the same. But so far this speculation has not been proved mathematically. In many cases, ε curves of physical quantities, such as the aforementioned Θ b − ε curve, can be applied to identify a suitable ε value to ensure that the series solutions Θ(ξ) converge in the space of the whole physical context. Therefore, it can be stated that the ε curve provides an effective way to study the influence of auxiliary parameter ε on the convergence region and rate for the series solutions.
By comparing the original Equation (47) with the auxiliary Equation (56), the following equation can be obtained: Substituting θ = π 2 into the Equation (81) leads to By solving the above equations, the optimal convergence parameter in the overall regions is derived as: When n = 1, the 1st-order approximation of θ b by IHAM is calculated by: For example, the 1st-order approximation of Θ b with ћ = −1, ε = 8/π 2 , is In the case of n = 1, ћ = −1, when ε = 0, ε = 1, and ε = 8/π 2 respectively, the 1st-order Θ  Table 2.    When the relative error ∆ ≤ 1%, the convergent region is: From Table 2, when ε = 0 the convergence region of Θ b is iteratively calculated to be 0 < α ≤ 0.3. When ε = 8/π 2 , the convergence region of Θ b is 0 < α ≤ 1.3. It can be stated that the improved homotopy analysis method expands the convergence region by 333% for this specific case.
In the case of n = 1, ћ = −1, the 1st-order solution Θ to rotation angle at the free end versus the point force by the improved homotopy analysis method, and the exact solution Θ b are compared in Figure 2, from which the optimal control parameter of the convergence region for the improved homotopy analysis method is identified to be ε = 8/π 2 . From Table 2, the series solutions to rotation angle converge to the exact solution and the relative error is less than 1% for 0 < α ≤ 0.3 when ε increases from 0 to 1. Thus, when 0 < α ≤ 0.3, the effective region of ε is specified to be R s = [0, 1].
In the case of n = 3, ε = 8/π 2 , the 3rd-order solution Θ to the rotation angle at the free end by IHAM vs. the point force curve is compared to the exact solutions Θ b in Figure 3. It can be stated that the improved homotopy analysis method has the optimal control parameters ћ for the convergence rate for an arbitrary nth-order approximation solution. For example, when n = 3, the optimal control parameter is identified to be ћ = −0.5.

Comparison between IHAM and HAM
The 30th-order approximation of θ b with n = 3, ћ = −0.5 by HAM [44] can be calculated as: When the relative error ∆ ≤ 1%, the convergent regions of the 30th-order approximation by HAM [44] are: The 30th-order approximation of θ b with ε = 8/π 2 , ћ = −2.25 by IHAM is When the relative error ∆ ≤ 1%, the convergent regions of the 30th-order approximation by IHAM are: Table 4 compares the 30th-order solutions by IHAM and HAM, from which it can be seen that the improved homotopy analysis method extends the convergence region from α ∈ [0, 1.3] by the homotopy analysis method to ћ = −0.5. Table 4. Thirtieth-order solutions Θ

•
Discrepancy between the linear solution and the exact solution is very large; • Thirtieth-order solution by IHAM with ε = 8/π 2 approaches the exact solution much closer than the 30th-order solution does by HAM with ε = 0; • The convergence region and rate have been substantially improved by homotopy-Páde approximation; • And again the Páde approximation solution Pade [15,15] by IHAM is closer to the exact solution than is the Páde approximation solution Pade [11,11] by HAM.

Discussion about the Relative Errors
In this section, the three iterations of α = 2.0 are taken as an example to discuss the effect of the control parameters ε and ћ on the relative errors.

Conclusions
The advance of the currently proposed IHAM from its prototype lies in the introduction of a nonlinear differential equation when constructing the homotopy equation. The solution by the homotopy nonlinear differential equation can approximate the exact solution of strong nonlinear differential equation at a higher convergence rate than the solution by the linear differential equation used in traditional homotopy analysis does.
As an application, the explicit solutions to the rotation angle for a large deformation of a cantilever beam under a point load at the free end are derived by IHAM. Current solutions with different parameters match well with the exact solutions from elliptical integrals. They are explicit to be simple and straightforward compared to the implicit exact solution from the elliptical integrals that are required to solve a transcendental equation. It can be stated that the current solution is easy to calculate with the explicit polynomial expressions, which facilitates practical engineering applications with low requirements on the calculation and computation.
By comparing the exact solution to the rotation angle of a large deformation of a cantilever beam, the effectiveness of the IHAM in solving strong nonlinear differential equations is demonstrated. The current improved homotopy analysis method has great superiority over the traditional homotopy analysis method in view of the rate and region of convergence, and the error of the solution by IHAM against the exact solution is reduced.
In another parallel study, the analytical expressions of vertical and horizontal displacements of the large deformation cantilever beam are deduced using the improved homotopy analysis method. It has been published in another paper: "EXPLICIT SOLU-TION TO LARGE DEFORMATION OF CANTILEVER BEAM BY IMPROVED HOMOTOPY ANALYSIS METHOD II: VERTICAL AND HORIZONTAL DISPLACEMENTS".