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

Abstract

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 firstly 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 significantly 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 deflection 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 significant 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.

1. 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 non-linear 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

$$\frac{\mathrm{d}x}{\mathrm{d}t}=A\left(t\right)x$$

(1)

where $A\left(t\right)$ is a time periodic matrix, Lyapunov [11] introduced an artificial parameter $m$. Thus, Equation (1) becomes

$$\frac{\mathrm{d}x}{\mathrm{d}t}=mA\left(t\right)x$$

(2)

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 $\delta$-expansion method by Karmishin et al. [12]. On the other hand, both the artificial parameter method and the $\delta$-expansion method fail to provide a principal to identify the domain and specify the artificial parameter or $\delta$. Similar to perturbation techniques, both the artificial small parameter and the $\delta$-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 $\mathrm{ћ}$ 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;
• It controls effectively the convergence of approximation series and adjusts convergence regions when necessary;
• It efficiently approximates a non-linear problem by choosing different sets of base functions.

As a proficient technique to solve strong nonlinear equations, HAM has been successfully employed to solve various types of nonlinear problems over the past two decades [23,24,25,26,27,28,29,30,31,32,33]. For example, the homotopy analysis method has been applied to solve nonlinear eigenvalue problems [19] and various nonlinear evolution equations [34,35,36,37,38,39].

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.

2. Improved Homotopy Analysis Method

2.1. 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. Improved homotopy analysis method vs. homotopy analysis method.

HAM	IHAM
$\mathrm{Original} \mathcal{N}\left(\xi \right)=0$ and $\mathrm{Auxiliary} \mathcal{L}\left(\xi \right)=0$ $\Rightarrow \mathrm{High} \mathrm{order} \mathcal{H}\left(\xi ;q\right)=0$	$① \mathrm{Original} \mathcal{N}\left(\xi \right)=0$ $\mathrm{and} \mathrm{Auxiliary} {\mathcal{N}}_0\left(\xi \right)=0$ $\Rightarrow \mathrm{New} \overline{\mathcal{N}}\left(\xi ;q\right)=0$
	$② \mathrm{New} \overline{\mathcal{N}}\left(\xi ;q\right)=0$ and $\mathrm{Auxiliary} \mathcal{L}\left(\xi \right)=0$ $\Rightarrow \mathrm{High} \mathrm{order} \overline{\mathcal{H}}\left(\xi ;q\right)=0$

Notes: $\mathcal{N}$—nonlinear operator, $\mathcal{L}$—lonear operator, ${\mathcal{N}}_0$—auxiliary nonlinear operator, $\overline{\mathcal{N}}$—new nonlinear operator, $\mathcal{H}$—homotopy operator, $\xi$—independent variable, $q$—embedded parameter.

lists the formulation steps of IHAM and HAM for comparison.

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

$$\mathcal{N}\left[w\left(\xi \right)\right]=0:w^{″}\left(\xi \right)+f\left[w\left(\xi \right)\right]=0,\xi \in \left[0,1\right]$$

(3)

Here $f\left(w\right)$ is a derivable function with the initial values of $f_w\left(0\right)=c_0$, $f_w^{\prime }\left(0\right)=c_1$, and $\frac{1}{2!}{f}_w^{″}\left(0\right)=c_2$. In compliance with the boundary conditions, we have

$$w\left(0\right)=0, w^{\prime }\left(1\right)=0$$

(4)

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.

Linearizing Equation (4) yields

$$\mathcal{L}\left[w\left(\xi \right)\right]=0:w^{″}\left(\xi \right)+c_0+c_{1}w\left(\xi \right)=0,\xi \in \left[0,1\right]:$$

(5)

with the concomitant boundary conditions

$$w\left(0\right)=0,w^{\prime }\left(1\right)=0$$

(6)

The auxiliary nonlinear operator is chosen as

$${\mathcal{N}}_0\left[w\left(\xi \right)\right]=0:w^{″}\left(\xi \right)+c_0+c_{1}w\left(\xi \right)+c_2\epsilon w^2\left(\xi \right)=0,\xi \in \left[0,1\right]:$$

(7)

wherein $\delta \in \left[0,1\right]$.

2.2. Traditional Homotopy Analysis Method

2.2.1. Zero-Order Deformation Equation

For comparison, nonlinear boundary value Equations (3) and (4) are solved firstly by HAM. When $w\left(\xi \right)$ is expanded into a power series of elementary functions in terms of $\xi$, we have

$$w\left(\xi \right)={\displaystyle \sum _{k=1}^{+\infty }{a}_k{\xi }^k}$$

(8)

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

$$w_0\left(\xi \right)=\frac{a_1}{2}\left(2-\xi \right)\xi$$

(9)

which will be specified as an initial guess solution for $w\left(\xi \right)$.

The homotopy analysis method is based on assumptive continuous mapping of $w\left(\xi \right)\to \mathsf{\Psi }\left(\xi ,q\right)$. When the embedded variable $q$ increases from 0 to 1, $\mathsf{\Psi }\left(\xi ,q\right)$ changes from the initial guess solution $w_0\left(\xi \right)$ to the exact solution $w\left(\xi \right)$. According to the Equations (3) and (8), the below auxiliary linear operator is selected

$$\mathrm{L}\left[\mathsf{\Psi }\left(\xi ,q\right)\right]=\frac{{\partial }^2\mathsf{\Psi }\left(\xi ,q\right)}{\partial {\xi }^2}$$

(10)

It satisfies

$$\mathcal{L}\left(C_1+C_2\xi \right)=0$$

(11)

where $C_1$ and $C_2$ are coefficients. From Equation (3), the following nonlinear operator is specified

$$\mathrm{N}\left[\mathsf{\Psi }\left(\xi ,q\right)\right]=\frac{{\mathrm{d}}^2\mathsf{\Psi }\left(\xi ,q\right)}{\mathrm{d}{\xi }^2}+f\left[q\mathsf{\Psi }\left(\xi ,q\right)\right]+c_1\left(1-q\right)\mathsf{\Psi }\left(\xi ,q\right)$$

(12)

Denoting a non-zero auxiliary parameter as $\mathrm{ћ}$, a more general form of homotopy is then constructed as

$$H\left[\mathsf{\Psi }\left(\xi ,q\right)\right]=\left(1-q\right)\mathrm{L}\left[\mathsf{\Psi }\left(\xi ,q\right)-w_0\left(\xi \right)\right]-q\mathrm{ћ}\mathrm{N}\left[\mathsf{\Psi }\left(\xi ,q\right)\right]$$

(13)

Setting $\left[\mathsf{\Psi }\left(\xi ,q\right)\right]$ = 0, a family of equations are derived as follows

$$\left(1-q\right)\mathrm{L}\left[\mathsf{\Psi }\left(\xi ,q\right)-w_0\left(\xi \right)\right]=q\mathrm{ћ}\mathrm{N}\left[\mathsf{\Psi }\left(\xi ,q\right)\right]$$

(14)

Equation (14) is subject to the boundary conditions

$$\mathsf{\Psi }\left(0,q\right)=0, {\mathsf{\Psi }}^{\prime }\left(1,q\right)=0$$

(15)

where $q\in \left[0,1\right]$ is an embedded parameter.

When $q=0$, Equations (14) and (15) yield

$$\mathsf{\Psi }\left(\xi ,0\right)=w_0\left(\xi \right)$$

(16)

When $q=1$, Equations (14) and (15) are exactly the same as Equations (3) and (4) provided that

$$\mathsf{\Psi }\left(\xi ,1\right)=w\left(\xi \right)$$

(17)

Therefore, as the embedded parameter $q$ increases from 0 to 1, $\mathsf{\Psi }\left(\xi ,q\right)$ evolves from the initial guess $w_0\left(\xi \right)$ to the solution $w\left(\xi \right)$. For brevity, Equations (14) and (15) are defined as zero-order deformation equations.

Assume that $\mathsf{\Psi }\left(\xi ,q\right)$ is analytic in the domain of $q\in \left[0,1\right]$, so that the deformation derivatives exist:

$$w_0^{\left[n\right]}\left(\xi \right)=\frac{\partial {\mathsf{\Psi }}^n\left(\xi ,q\right)}{\partial q^n}{|}_{q=0}$$

(18)

$\mathsf{\Psi }\left(\xi ,q\right)$ can then be expanded in the Maclaurin series of $q$ as follows:

$$\mathsf{\Psi }\left(\xi ,q\right)=w_0\left(\xi \right)+\sum _{n=1}^{+\infty }{w}_n\left(\xi \right)q^n$$

(19)

where

$$w_n\left(\xi \right)={\left.\frac{1}{n!}\frac{{\partial }^n\mathsf{\Psi }\left(\xi ,q\right)}{\partial q^n}\right|}_{q=0}$$

(20)

It is worth mentioning that the auxiliary parameter $\mathrm{ћ}$ in Equation (14) determines the convergence regions. Assuming that $\mathrm{ћ}$ is properly chosen, so that all of these Maclaurin series are convergent at $q=1$. Thus, at $q=1$ Equation (17) gives

$$w\left(\xi \right)=w_0\left(\xi \right)+\sum _{n=1}^{+\infty }{w}_n\left(\xi \right)$$

(21)

Equation (21) bridges the initial guess solution $w_0\left(\xi \right)$ to the exact solution $w\left(\xi \right)$. The result for the $n$th-order approximation is given by:

$$w\left(\xi \right)\approx w_0\left(\xi \right)+{\sum }_{m=1}^n{w}_m\left(\xi \right)$$

(22)

2.2.2. High-Order Deformation Equation

Define a vector

$$w_n=\left\{w_0\left(\xi \right),w_1\left(\xi \right),w_2\left(\xi \right),\dots , w_n\left(\xi \right)\right\}$$

(23)

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\left(\xi \right)$,

$$\mathrm{L}\left[w_n\left(\xi \right)-{\chi }_n{w}_{n-1}\left(\xi \right)\right]=\mathrm{ћ}{R}_n\left(\xi ,q\right)$$

(24)

in conjunction with the boundary conditions

$$w_n\left(0\right)=0,w_n^{\prime }\left(1\right)=0$$

(25)

where

$$R_n={\left.\frac{1}{\left(n-1\right)!}\frac{{\partial }^{n-1}\mathrm{N}\left[\mathsf{\Psi }\left(\xi ,q\right),q\right]}{\partial q^{n-1}}\right|}_{q=0}$$

(26)

and

$${\chi }_n=\left\{\begin{array}{c}0,\hspace{1em}n\le 1\hfill \\ 1,\hspace{1em}n>1\hfill \end{array}\right.$$

(27)

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.

2.3. Improved Homotopy Analysis Method

2.3.1. 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:

$$\begin{array}{c}\overline{\mathrm{N}}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]=\frac{{\mathrm{d}}^2\overline{\mathsf{\Psi }}\left(\xi ,q\right)}{\mathrm{d}{\xi }^2}+f\left[q\cdot \overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]+c_1\left(1-q\right)\overline{\mathsf{\Psi }}\left(\xi ,q\right)\\ +c_2\epsilon \left(1-q^2\right){\overline{\mathsf{\Psi }}}^2\left(\xi ,q\right)\end{array}$$

(28)

The improved homotopy analysis method replaces the nonlinear operator $\mathcal{N}$ in Equation (12) with the new nonlinear operator $\overline{\mathcal{N}}$ in Equation (28).

When $\epsilon =0$, the new nonlinear operator $\overline{\mathcal{N}}$ in Equation (28) degenerates to the nonlinear operator $\mathcal{N}$ in Equation (12).

When $q=1$, the new nonlinear differential equations $\overline{\mathcal{N}}=0$ becomes the original differential Equation (3) for $\mathcal{N}=0$.

When $q=0$, the nonlinear differential equations $\overline{\mathcal{N}}=0$ is the simple differential Equation (7) for $\mathcal{N}=0$.

2.3.2. Construction of High-Order Homotopy Equation

With the new nonlinear operator $\overline{\mathcal{N}}$ in Equation (28) and the selected linear operator $\mathcal{L}$ in Equation (10), the homotopy operator is defined as follows,

$$\overline{\mathrm{H}}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]=\left(1-q\right)\mathrm{L}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)-w_0\left(\xi \right)\right]-q\mathrm{ћ}\overline{\mathrm{N}}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]$$

(29)

The improved homotopy analysis method replaces the homotopy operator $\mathcal{H}$ in Equation (13) with the new homotopy operator $\overline{\mathcal{H}}$ in Equation (29).

When $\epsilon =0$, the new homotopy operator $\overline{\mathcal{H}}$ in Equation (29) degenerates to the homotopy operator $\mathcal{H}$ in Equation (13).

With definition of a vector

$${\overline{w}}_n=\left\{w_0\left(\xi \right),{\overline{w}}_1\left(\xi \right),{\overline{w}}_2\left(\xi \right),\cdots ,{\overline{w}}_n\left(\xi \right)\right\}$$

(30)

a new homotopy equation of high order is similarly constructed

$$\mathrm{L}\left[{\overline{w}}_n\left(\xi \right)-{\chi }_n{\overline{w}}_{n-1}\left(\xi \right)\right]=\mathrm{ћ}{\overline{R}}_n\left({\overline{w}}_{n-1},\xi \right)$$

(31)

in conjunction with the boundary conditions

$$w\left(0\right)=0,w^{\prime }\left(1\right)=0$$

(32)

where

$${\overline{R}}_n\left({\overline{w}}_{n-1},\xi \right)={\left.\frac{1}{\left(n-1\right)!}\frac{{\partial }^{n-1}\overline{\mathrm{N}}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]}{\partial q^{n-1}}\right|}_{q=0}$$

(33)

and

$${\chi }_n=\left\{\begin{array}{c}0,\hspace{1em}n\le 1\hfill \\ 1,\hspace{1em}n>1\hfill \end{array}\right.$$

(34)

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 $n$th-order approximation are given by

$$\overline{w}\left(\xi \right)\approx w_0\left(\xi \right)+\sum _{m=1}^n{\overline{w}}_m\left(\xi \right)$$

(35)

2.3.3. Convergence Theorem

Theorem 1 (Convergence Theorem).

If the series

$$w_0\left(\xi \right)+\sum _{n=1}^{+\infty }{\overline{w}}_n\left(\xi \right)$$

(36)

converges, where${\overline{w}}_n\left(\xi \right)$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

$$\overline{S}\left(\xi \right)=w_0\left(\xi \right)+\sum _{n=1}^{+\infty }{\overline{w}}_n\left(\xi \right)$$

(37)

being the convergent series, a necessary condition for the series to converge is:

$$\underset{n\to +\infty }{lim}{\overline{w}}_n\left(\xi \right)=0$$

(38)

From Equations (31), (34) and (38), we have

$$\begin{gathered} \mathrm{ћ}\sum _{n=1}^{+\infty }{\overline{R}}_n\left({\overline{w}}_{n-1},\xi \right) \\ =\sum _{n=1}^{+\infty }\mathcal{L}\left[{\overline{w}}_n\left(\xi \right)-{\chi }_n{\overline{w}}_{n-1}\left(\xi \right)\right] \\ =\mathcal{L}\left\{\sum _{n=1}^{+\infty }\left[{\overline{w}}_n\left(\xi \right)-{\chi }_n{\overline{w}}_{n-1}\left(\xi \right)\right]\right\} \\ =\mathcal{L}\left[{\overline{w}}_1+\left({\overline{w}}_2-{\overline{w}}_1\right)+\left({\overline{w}}_3-{\overline{w}}_2\right)+\dots +\left({\overline{w}}_n-{\overline{w}}_{n-1}\right)+\dots \right] \\ =\mathcal{L}\left[\underset{n\to +\infty }{lim}{\overline{w}}_n\left(\xi \right)\right]=0 \end{gathered}$$

Because $\mathrm{ћ}\ne 0$, so

$$\sum _{n=1}^{+\infty }{\overline{R}}_n\left({\overline{w}}_{n-1},\xi \right)=0$$

(39)

From Equation (33)

$$\sum _{n=1}^{+\infty }{\overline{R}}_n\left({\overline{w}}_{n-1},\xi \right)=\sum _{n=1}^{+\infty }{\left.\frac{1}{\left(n-1\right)!}\frac{{\partial }^{n-1}\overline{N}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]}{\partial q^{n-1}}\right|}_{q-0}=0$$

(40)

Generally, $\overline{\mathsf{\Psi }}\left(\xi ,q\right)$ does not satisfy the original nonlinear Equation (3). Letting

$$\overline{\delta }\left(\xi ,q\right)=\overline{\mathrm{N}}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]$$

(41)

represents the residual error of Equation (3), then

$$\overline{\delta }\left(\xi ,q\right)=0$$

corresponds to the exact solution of Equation (3). According to the above definition, the McLaughlin series of residual error $\overline{\delta }\left(\xi ;q\right)$ with respect to the embedded variable $q$ is:

$$\sum _{n=0}^{+\infty }{\left.\frac{q^n}{n!}\frac{{\partial }^n\overline{\delta }\left(\xi ,q\right)}{\partial q^n}\right|}_{q=0}=\sum _{n=0}^{+\infty }{\left.\frac{q^n}{n!}\frac{{\partial }^n\overline{\mathrm{N}}\left[\overline{\mathsf{\Psi }}\left(\xi ,q\right)\right]}{\partial q^n}\right|}_{q=0}$$

(42)

According to Equation (40), when $q$ = 1, Equation (42) becomes:

$$\overline{\delta }\left(\xi ,1\right)=\sum _{n=0}^{+\infty }{\left.\frac{1}{n!}\frac{{\partial }^n\overline{\delta }\left(\xi ,q\right)}{\partial q^n}\right|}_{q-0}=0$$

(43)

According to the definition of $\overline{\delta }\left(\xi ;q\right)$, when $q=1$, Equation (43) is the exact solution to the original Equation (3). Therefore, as long as the below series:

$$\mid w_0\left(\xi \right)+\sum _{n=1}^{+\infty }{\overline{w}}_n\left(\xi \right)$$

converges, there exists a solution to the original Equation (3).

Theorem 2.

If the series in Equation (36) converges, where${\overline{w}}_n\left(\xi \right)$satisfies the high-order deformation Equations (31) and (32), and the definitions in Equations (33) and (34) are true, there is:

$$\sum _{n=1}^{+\infty }{\overline{R}}_n\left({\overline{w}}_{n-1},\xi \right)=0$$

(44)

This theorem is self-explanatory with reference to Equation (40).

According to Theorems 1 and 2, only the initial guessing solution $w_0\left(\xi \right)$, auxiliary parameter $\mathrm{ћ}$, auxiliary linear operator $\mathcal{L}$, and auxiliary nonlinear operator $\mathcal{N}$ to ensure the convergence of series in Equation (36) should be derived.

3. Explicit Solution to Rotation Angle for Large Deformation of Cantilever Beams by IHAM

3.1. 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]

$$\frac{\mathrm{d}\theta }{\mathrm{d}s}=\frac{F}{EI}\left(l-x\right)$$

(45)

It complies to the boundary conditions:

$$\theta \left(0\right)=0,\hspace{1em}{\theta }^{\prime }\left(L\right)=0,$$

(46)

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, $\theta$ 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 $\xi =S/L$, the original Equation (45) becomes

$${\theta }^{″}+\alpha \cos \theta =0$$

(47)

in compliance to boundary conditions of

$$\theta \left(0\right)=0,\hspace{1em}{\theta }^{\prime }\left(1\right)=0,$$

(48)

The prime in Equations (47) and (48) denotes differentiation with respect to $\xi$, and

$$\alpha =\frac{F{L}^2}{EI}$$

(49)

The rotation angle of cross-section plane at the free end is denoted by ${\theta }_b=\theta \left(1\right)$.

For infinitesimal deformation, the below linear equation will be sufficient

$${\theta }^{″}+\alpha =0$$

(50)

which is subjected to the boundary conditions

$$\theta \left(0\right)=0,\hspace{1em}{\theta }^{\prime }\left(1\right)=0,$$

(51)

The corresponding solution is:

$$\theta \left(\xi \right)=\frac{\alpha }{2}\left(2-\xi \right)\xi$$

(52)

which gives the linear result of ration angle at the free end for $\xi =1$

$${\theta }_b=\frac{\alpha }{2}$$

(53)

3.2. Zero Order Deformation Equation

The nonlinear boundary value Equations (47) and (48) are then solved by IHAM to derive a explicit formula for ${\theta }_b$. To start with, $\theta \left(\xi \right)$ is expanded into the power series of $\xi$:

$$\theta \left(\xi \right)=\sum _{k=1}^{+\infty }{a}_k{\xi }^k$$

(54)

Equation (52) is then identified as an initial guess. Accordingly, the initial guess of $\theta \left(\xi \right)$ is expressed as:

$${\theta }_0\left(\xi \right)=\frac{\alpha }{2}\left(2-\xi \right)\xi$$

(55)

For simplicity, the below simplest auxiliary nonlinear differential equation is chosen

$${\theta }^{″}+\alpha \left(1-\frac{\epsilon }{2!}{\theta }^2\right)=0$$

(56)

which is subjected to the boundary conditions

$$\theta \left(0\right)=0,\hspace{1em}{\theta }^{\prime }\left(1\right)=0,$$

(57)

where $\epsilon$ is the correction coefficient of the convergent region, and $\epsilon \in \left[0,1\right]$.

When $\epsilon =0$, Equation (56) becomes the linear Equation (50). When $\epsilon =1$, Equation (56) becomes:

$${\theta }^{″}+\alpha \left(1-\frac{1}{2!}{\theta }^2\right)=0$$

(58)

In fact, the approximating second-order Taylor expansion for the cosine function in Equation (47) can be written as:

$$\cos \theta \approx 1-\frac{1}{2!}{\theta }^2$$

According to the original Equation (47) and the equivalent Equation (56), we define one nonlinear operator as:

$$\left[\psi \left(\xi ,q\right),q,\epsilon \right]=\frac{{\mathrm{d}}^2\psi \left(\xi ,q\right)}{\mathrm{d}{\xi }^2}+\alpha \cos \left[q\cdot \psi \left(\xi ,q\right)\right]-\frac{\alpha \epsilon }{2!}\left(1-q^2\right){[\psi \left(\xi ,q\right)]}^2$$

(59)

where $q\in \left[0,1\right]$ is an embedded parameter, $\psi \left(\xi ,q\right)$ is a function dependent on ξ and $q$, $\epsilon$ is a control parameter for the convergent region, and $\epsilon \in \left[0,1\right]$.

When $q=0$, there is:

$$\mathrm{N}\left[\psi \left(\xi ,0\right),0,\epsilon \right]=\frac{{\mathrm{d}}^2\psi \left(\xi ,0\right)}{\mathrm{d}{\xi }^2}+\alpha -\frac{\alpha \epsilon }{2!}{[\psi \left(\xi ,0\right)]}^2$$

If $q=1$, we have

$$\mathrm{N}\left[\psi \left(\xi ,1\right),1,\epsilon \right]=\frac{{\mathrm{d}}^2\psi \left(\xi ,1\right)}{\mathrm{d}{\xi }^2}+\alpha \cos \left[\psi \left(\xi ,1\right)\right]$$

Letting

$$\mathrm{L}\left[\psi \left(\xi ,q\right)\right]=\frac{{\partial }^2\psi \left(\xi ,q\right)}{\partial {\xi }^2}$$

(60)

Being an auxiliary linear operator, and $\mathrm{ћ}\in \left[-1,0\right)$ denotes a nonzero auxiliary parameter (convergence–control parameter), a homotopy is then constructed as:

$$\mathrm{H}\left[\psi \left(\xi ;q\right),q\right]=\left(1-q\right)\mathrm{L}\left[\psi \left(\xi ,q\right)-{\theta }_0\left(\xi \right)\right]-q\mathrm{hN}\left[\psi \left(\xi ,q\right),q,\epsilon \right]$$

(61)

When $q=0$, we have

$$\mathrm{H}\left[\psi \left(\xi ,0\right),0\right]=\mathrm{L}\left[\psi \left(\xi ,q\right)-{\theta }_0\left(\xi \right)\right]$$

If $q=1$, Equation (61) becomes:

$$\mathrm{H}\left[\psi \left(\xi ,1\right),1\left]=-\mathrm{hN}\right[\psi \left(\xi ,1\right),1,\epsilon \right]=-\mathrm{h}\left\{\frac{{\mathrm{d}}^2\psi \left(\xi ,1\right)}{\mathrm{d}{\xi }^2}+\alpha \cos \left[\psi \left(\xi ,1\right)\right]\right\}$$

Thus, by enforcing

$$\psi \left(0,q\right)=0,{\psi }_{\xi }^{\prime }\left(1,q\right)=0$$

We have a family of equations:

$$\left(1-q\right)\mathcal{L}\left[\psi \left(\xi ;q\right)-{\theta }_0\left(\xi \right)\right]=q\mathrm{ћ}𝓇\left[\psi \left(\xi ;q\right),q,\epsilon \right]$$

(62)

which are subjected to the boundary conditions of:

$$\psi \left(0,q\right)=0,{\psi }_{\xi }^{\prime }\left(1,q\right)=0$$

(63)

The prime in Equation (63) denotes differentiation with respect to $\xi$.

When $q=0$, Equation (62) incurs:

$$\psi \left(\xi ,0\right)={\theta }_0\left(\xi \right)$$

(64)

When $q=1$, Equations (62) and (63) are equivalent to the original Equations (47) and (48), provided that:

$$\psi \left(\xi ,1\right)=\theta \left(\xi \right)$$

(65)

Thus, as $q$ increases from 0 to 1, $\psi \left(\xi ,q\right)$ evolves from the known initial guess ${\theta }_0\left(\xi \right)$ to the solution $\theta \left(\xi \right)$ in Equations (47) and (48).

By expanding $\psi \left(\xi ;q\right)$ into a Taylor series of the embedded parameter $q$ and using Equation (64), we have:

$$\psi \left(\xi ;q\right)={\theta }_0\left(\xi \right)+\sum _{n=1}^{+\infty }{\theta }_n\left(\xi \right)q^n$$

(66)

where

$${\theta }_n={\left.\frac{1}{n!}\frac{{\partial }^n\psi \left(\xi ;q\right)}{\partial q^n}\right|}_{q=0}$$

Assuming that $\mathrm{ћ}$ and $\epsilon$ are properly chosen so that the series Equation (66) are convergent at $q=1$, according to Equation (65), we have the series:

$$\theta \left(\xi \right)={\theta }_0\left(\xi \right)+\sum _{n=1}^{+\infty }{\theta }_n\left(\xi \right)$$

(67)

And the governing equations of ${\theta }_n\left(\xi \right)$ can then be deduced from the zero-order deformation Equations (62) and (63).

3.3. 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 $n$th-order deformation equations as follows

$$\mathcal{L}\left[{\theta }_n\left(\xi \right)-{\chi }_n{\theta }_{n-1}\left(\xi \right)\right]=\mathrm{ћ}{R}_n\left({\theta }_0,{\theta }_1,\cdots ,{\theta }_{n-1}\right)$$

(68)

Equation (68) is subjected to the boundary conditions:

$${\theta }_{\pi }\left(0\right)=0,{\theta }_{\pi }^{\prime }\left(1\right)=0,$$

(69)

where

$$R_n={\left.\frac{1}{\left(n-1\right)!}\frac{{\partial }^{n-1}\mathcal{N}\left[\psi \left(\xi ;q\right),q,\epsilon \right]}{\partial q^{n-1}}\right|}_{q=0}$$

(70)

and

$${\chi }_n=\left\{\begin{array}{c}0,\hspace{1em}n\le 1\hfill \\ 1,\hspace{1em}n>1\hfill \end{array}\right.$$

(71)

Combining Equations (59) and (70), we have:

$$R_1={\theta }_0^{″}\left(\xi \right)+\alpha -\frac{1}{2}\epsilon \alpha {\theta }_0^2\left(\xi \right)$$

(72a)

$$R_2={\theta }_1{}^{″}\left(\xi \right)-\alpha \epsilon {\theta }_0\left(\xi \right){\theta }_1\left(\xi \right)$$

(72b)

$$R_3={\theta }_2^{″}\left(\xi \right)-\frac{1}{2}\alpha {\theta }_0^2\left(\xi \right)+\frac{1}{2}\alpha \epsilon {\theta }_0^2\left(\xi \right)-\frac{1}{2}\alpha \epsilon {\theta }_1{}^2\left(\xi \right)-\alpha \epsilon {\theta }_0\left(\xi \right){\theta }_2\left(\xi \right)$$

(72c)

$$R_4={\theta }_3{}^{″}\left(\xi \right)-\alpha {\theta }_0\left(\xi \right){\theta }_1\left(\xi \right)+\alpha \epsilon {\theta }_0\left(\xi \right){\theta }_1\left(\xi \right)-\alpha \epsilon {\theta }_1\left(\xi \right){\theta }_2\left(\xi \right)-\alpha \epsilon {\theta }_0\left(\xi \right){\theta }_3\left(\xi \right)$$

(72d)

$$\begin{array}{l}{R}_5={\theta }_4^{″}\left(\xi \right)-\alpha {\theta }_0\left(\xi \right){\theta }_2\left(\xi \right)-\frac{1}{2}\alpha {\theta }_1{}^2\left(\xi \right)+\frac{\alpha }{24}{\theta }_0^4\left(\xi \right)+\alpha \epsilon {\theta }_0\left(\xi \right){\theta }_2\left(\xi \right)\\ -\alpha \epsilon {\theta }_1\left(\xi \right){\theta }_3\left(\xi \right)-\alpha \epsilon {\theta }_0\left(\xi \right){\theta }_4\left(\xi \right)+\frac{1}{2}\alpha \epsilon {\theta }_1{}^2\left(\xi \right)-\frac{1}{2}\alpha \epsilon {\theta }_2^2\left(\xi \right)\end{array}$$

(72e)

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 $\epsilon$ about nonlinear operator ${\mathcal{N}}_0$ and the auxiliary parameter $\mathrm{ћ}$ about linear operator $\mathcal{L}$, which control the convergence region and the rate of the IHAM solution series.

4. Calculation Results

For simplicity, the rotation angle with dimension one is defined as:

$${\mathsf{\Theta }}_b=\frac{{\theta }_b}{\pi /2}$$

(73)

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 ${\theta }_0\left(\xi \right)$ and auxiliary linear operator $\mathcal{L}$. It would be helpful to put forward the so-called solution expression principle: the initial guess solution ${\theta }_0\left(\xi \right)$ and the solution ${\theta }_1\left(\xi \right),{\theta }_2\left(\xi \right), \dots$ of the high-order deformation equation should not violate the solution expression.

Provided that the initial guess solution ${\theta }_0\left(\xi \right)$ and the auxiliary linear operator $\mathcal{L}$ are chosen, the homotopy analysis method still affords flexibility in identifying the value of the auxiliary parameter $\mathrm{ћ}$. Superior to all the traditional analytical methods, HAM always offers a series of solutions with auxiliary parameter $\mathrm{ћ}$. The auxiliary parameter $\mathrm{ћ}$ dominates the convergence region and rate of the series of solutions. By selecting an appropriate $\mathrm{ћ}$ value, the convergence region and the rate of the series of solutions can be effectively controlled.

Similar to the auxiliary parameter $\mathrm{ћ}$, even if the initial guessing solution ${\theta }_0\left(\xi \right)$, auxiliary linear operator $\mathcal{L}$ and auxiliary parameter $\mathrm{ћ}$ are specified, again IHAM is flexible with the value of auxiliary parameter $\epsilon$. Different from the original HAM, the IHAM gets a series of solutions with the variable auxiliary parameter $\epsilon$, which also affects the convergence region and rate of the series of solutions. By selecting the appropriate $\epsilon$ 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 $\mathrm{ћ}$ and $\epsilon$.

4.1. Exact Solution to Rotation Angle

From Equation (47), the transcendental equation solved ${\theta }_b$ is as follows

$$\sqrt{\alpha }=K\left(\mu \right)-F\left(\varphi ,\mu \right)$$

(74)

where

$$F\left(\varphi ,\mu \right)={\int }_0^{\varphi }\frac{\mathrm{d}t}{\sqrt{1-{\mu }^2{\sin }^{2}t}}$$

(75a)

$$K\left(\mu \right)={\int }_0^{\pi /2}\frac{\mathrm{d}t}{\sqrt{1-{\mu }^2{\sin }^{2}t}}$$

(75b)

$F\left(\varphi ,\mu \right)$ is the first type of elliptic integral and $K\left(\mu \right)$ is the first complete elliptic integral. Therein

$$\varphi =\arcsin \frac{1}{\sqrt{2}\mu }$$

(76)

$$\mu =\sqrt{\frac{1+\sin {\theta }_b}{2}}$$

(77)

According to Equations (53) and (73), the linear solution to the rotation angle at the free end with dimension 1 is:

$${\mathsf{\Theta }}_b=\frac{\alpha }{\pi }$$

(78)

It can be verified that the convergent region of the linear solution (78) is as follows if the relative error $\mathsf{\Delta }\le 1\%$ is required.

$$\alpha \in \left[0,0.33\right],\hspace{1em}{\mathsf{\Theta }}_b\in \left[0,0.11\right],$$

(79)

4.2. Effect of Auxiliary Parameter $\epsilon$ on Convergence

It is expected that the improved homotopy analysis method derives a family of series solutions with the variable auxiliary parameter $\epsilon$. The efficiency of the IHAM lies in the appropriate identification for the value of $\epsilon$ 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 $\epsilon$. Therefore, the curve of these physical quantities can be drawn with respect to $\epsilon$ by regarding the auxiliary parameter $\epsilon$ as an independent variable. For example, below ${\mathsf{\Theta }}_b$ is a physical quantity

$${\mathsf{\Theta }}_b={\left.\mathsf{\Theta }\left(\xi ,\epsilon \right)\right|}_{\xi -1}$$

(80)

If ${\mathsf{\Theta }}_b$ is a function of, $\epsilon , {\mathsf{\Theta }}_b-\epsilon$ can be identified. According to Theorem 1 or Theorem 2, all convergent series of ${\mathsf{\Theta }}_b$ given by different values of $\epsilon$ 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 ${\mathsf{\Theta }}_b-\epsilon$ curve to imply the converge region of $\epsilon$ represented by $R_{\epsilon }$. For simplicity, such a curve is called a $\epsilon$ curve, and the corresponding region $R_{\epsilon }$ is the effective region of $\epsilon$. Obviously, if the value of $\epsilon$ falls in the valid region, the series solution will converge without fail. Even though ${\mathsf{\Theta }}_b$ does not make any physical sense, the corresponding $\epsilon$ curve can still be attempted. The more $\epsilon$ curves that are drawn, the more specific value of $\epsilon$ will be identified. If the initial guess solution ${\theta }_0\left(\xi \right)$, auxiliary linear operator $\mathcal{L}$ and the auxiliary parameter $\mathrm{ћ}$ are given, the effective region of $\epsilon$ for different physical quantities is usually almost the same. But so far this speculation has not been proved mathematically. In many cases, $\epsilon$ curves of physical quantities, such as the aforementioned ${\mathsf{\Theta }}_b-\epsilon$ curve, can be applied to identify a suitable $\epsilon$ value to ensure that the series solutions $\mathsf{\Theta }\left(\xi \right)$ converge in the space of the whole physical context. Therefore, it can be stated that the $\epsilon$ curve provides an effective way to study the influence of auxiliary parameter $\epsilon$ 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:

$$\cos \theta =1-\frac{\epsilon }{2!}{\theta }^2$$

(81)

Substituting $\theta =\frac{\pi }{2}$ into the Equation (81) leads to

$$0=1-\frac{\epsilon }{2!}{\left(\frac{\pi }{2}\right)}^2$$

(82)

By solving the above equations, the optimal convergence parameter in the overall regions is derived as:

$$\epsilon =\frac{8}{{\pi }^2}$$

(83)

When $n=1$, the 1st-order approximation of ${\theta }_b$ by IHAM is calculated by:

$${\theta }_b=0.5\alpha +0.045833\mathrm{ћ}\epsilon {\alpha }^3$$

(84)

For example, the 1st-order approximation of ${\mathsf{\Theta }}_b$ with $\mathrm{ћ}=-1,\epsilon =8/{\pi }^2$, is

$${\mathsf{\Theta }}_b=0.31831\alpha -0.023651{\alpha }^3$$

(85)

In the case of $n=1, \mathrm{ћ}=-1$, when $\epsilon =0, \epsilon =1$, and $\epsilon =8/{\pi }^2$ respectively, the 1st-order ${\mathsf{\Theta }}_b^{\left[n=1\right]}$ by IHAM and the exact solution ${\mathsf{\Theta }}_b$ are compared in

Table 2. First-order numerical solutions ${\mathsf{\Theta }}_b^{\left[n=1\right]}$ by IHAM vs. exact solutions ${\mathsf{\Theta }}_b$.

$\mathit{\alpha }$	Exact Solution	$\mathit{\epsilon }=0$ (HAM [48])		$\mathit{\epsilon }=8/{\mathit{\pi }}^2$ (IHAM)		$\mathit{\epsilon }=0$ (IHAM)
	${\mathsf{\Theta }}_{\mathit{b}}$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$
0	0	0	0	0	0	0	0
0.1	0.031802	0.031831	0.0916%	0.031807	0.0172%	0.031802	0.000178%
0.2	0.063430	0.063662	0.365%	0.063473	0.0669%	0.063429	0.00284%
0.3	0.094719	0.095493	0.817%	0.094854	0.143%	0.094705	0.0143%
0.4	0.12551	0.12732	1.44%	0.12581	0.237%	0.12546	0.0446%
0.5	0.15567	0.15915	2.24%	0.15620	0.336%	0.15551	0.107%
0.6	0.18509	0.19099	3.19%	0.18588	0.426%	0.18468	0.219%
0.7	0.21366	0.22282	4.29%	0.21470	0.488%	0.21281	0.399%
0.8	0.24132	0.25465	5.52%	0.24254	0.505%	0.23971	0.667%
0.9	0.26801	0.28648	6.89%	0.26924	0.458%	0.26521	1.05%
1.0	0.29371	0.31831	8.38%	0.29466	0.324%	0.28913	1.56%
1.1	0.31839	0.35014	9.97%	0.31866	0.0857%	0.31130	2.22%
1.2	0.34206	0.38197	11.7%	0.34110	0.279%	0.33155	3.07%
1.3	0.36473	0.41380	13.5%	0.36184	0.791%	0.34970	4.12%
1.4	0.38641	0.44563	15.3%	0.38074	1.47%	0.36557	5.39%
1.5	0.40714	0.47746	17.3%	0.39764	2.33%	0.37899	6.92%
1.6	0.42685	0.50930	19.2%	0.41242	3.40%	0.38978	8.71%
1.7	0.44587	0.54113	21.4%	0.42493	4.70%	0.39777	10.8%
1.8	0.46394	0.57296	23.5%	0.43502	6.23%	0.40279	13.2%
1.9	0.48120	0.60479	25.7%	0.44257	8.03%	0.40465	15.9%
2	0.49768	0.63662	27.9%	0.44741	10.1%	0.40319	19.0%
$\infty$	1	1	0	1	0	1	0

.

When the relative error $\mathsf{\Delta }\le 1\%$, the convergent region is:

$$\alpha \in \left[0,1.3\right],{\hspace{1em}\mathsf{\Theta }}_b\in \left[0,0.36\right]$$

(86)

From Table 2, when $\epsilon =0$ the convergence region of ${\mathsf{\Theta }}_b$ is iteratively calculated to be $0<\alpha \le 0.3$. When $\epsilon =8/{\pi }^2$, the convergence region of ${\mathsf{\Theta }}_b$ is $0<\alpha \le 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,\mathrm{ћ}=-1$, the 1st-order solution ${\mathsf{\Theta }}_b^{\left[n-1\right]}$ to rotation angle at the free end versus the point force by the improved homotopy analysis method, and the exact solution ${\mathsf{\Theta }}_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 $\epsilon =8/{\pi }^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<\alpha \le 0.3$ when $\epsilon$ increases from 0 to 1. Thus, when $0<\alpha \le 0.3$, the effective region of $\epsilon$ is specified to be $R_s=\left[0,1\right]$.

4.3. Effect of Auxiliary Parameter $\mathrm{ћ}$ on Convergence

When $n=3,\epsilon =8/{\pi }^2$, the 3rd-order approximation of ${\theta }_b$ by IHAM is:

$$\begin{array}{c}{\theta }_b=0.5\alpha +\left(0.12014\mathrm{ћ}+0.11145{\mathrm{ћ}}^2+0.037151{\mathrm{ћ}}^3\right){\alpha }^3\hfill \\ +\left(0.016031{\mathrm{ћ}}^2+0.010687{\mathrm{ћ}}^3\right){\alpha }^5+9.5862\times {10}^{-4}{\mathrm{ћ}}^3{\alpha }^7\hfill \end{array}$$

(87)

From Equation (87), the 3rd-order approximation of ${\mathsf{\Theta }}_b$ with $\mathrm{ћ}=-1/2$ is calculated as:

$${\mathsf{\Theta }}_b=0.31831\alpha -2.3458\times {10}^{-2}{\alpha }^3+1.7010\times {10}^{-3}{\alpha }^5-7.6284\times {10}^{-5}{\alpha }^7$$

(88)

When the relative error $\mathsf{\Delta }\le 1\%$, the convergent region of the 3rd-order approximation by IHAM is:

$$\alpha \in \left[0,2.0\right],\hspace{1em}{\mathsf{\Theta }}_b\in \left[0,0.49\right],$$

(89)

In the case of $n=3$, $\epsilon =8/{\pi }^2$, when $\mathrm{ћ}=-0.6$, $\mathrm{ћ}=-0.4$ and $\mathrm{ћ}=-0.5$ respectively, the 3rd-order ${\mathsf{\Theta }}_b^{\left[n=3\right]}$ by IHAM are compared and with the exact solution in

Table 3. Third-order solutions ${\mathsf{\Theta }}_b^{\left[n=3\right]}$ by IHAM with different $\mathrm{ћ}$ and fixed $\epsilon =8/{\pi }^2$.

$\mathit{\alpha }$	Exact Solution	$\mathbf{ћ}=-0.6$		$\mathbf{ћ}=-0.5$		$\mathbf{ћ}=-0.4$
	${\mathsf{\Theta }}_{\mathit{b}}$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$
0	0	0	0	0	0	0	0
0.2	0.063430	0.063459	0.0452%	0.063475	0.0702%	0.063496	0.104%
0.4	0.12551	0.12572	0.163%	0.12584	0.261%	0.12601	0.395%
0.6	0.18509	0.18566	0.306%	0.18605	0.518%	0.18660	0.814%
0.8	0.24132	0.24231	0.411%	0.24318	0.771%	0.24441	1.28%
1.0	0.29371	0.29493	0.416%	0.29648	0.943%	0.29871	1.71%
1.2	0.34206	0.34300	0.275%	0.34539	0.975%	0.34895	2.01%
1.4	0.38641	0.38626	0.0408%	0.38961	0.827%	0.39472	2.14%
1.6	0.42685	0.42461	0.547%	0.42900	0.479%	0.43580	2.07%
1.8	0.46394	0.45810	1.26%	0.46362	0.0691%	0.47216	1.77%
2.0	0.49768	0.48666	2.21%	0.49362	0.815%	0.50391	1.25%
2.2	0.52847	0.50998	3.50%	0.51913	1.77%	0.53127	0.530%
2.4	0.55661	0.52714	5.29%	0.54011	2.96%	0.55449	0.382%
2.6	0.58237	0.53627	7.92%	0.55613	4.51%	0.57375	1.48%
2.8	0.60600	0.53404	11.9%	0.56612	6.58%	0.58908	2.79%
3	0.62772	0.51507	17.9%	0.56805	9.50%	0.60016	4.39%
$\infty$	1	1	0	1	0	1	0

.

By IHAM with only three calculation iterations, the convergence region of ${\mathsf{\Theta }}_b$ is derived to be $0<\alpha \le 0.6$, $0<\alpha \le 2.0$, and $0<\alpha \le 1.6$ when $\mathrm{ћ}=-0.4$, $\mathrm{ћ}=-0.5$, and $\mathrm{ћ}=-0.6$, respectively. It is observed from Table 3 that the convergence region expands by 233% when $\mathrm{ћ}=-0.4$ changes to $\mathrm{ћ}=-0.5$, whereas it expands by 25% when $\mathrm{ћ}=-0.6$ changes to $\mathrm{ћ}=-0.5$.

In the case of $n=3$, $\epsilon =8/{\pi }^2$, the 3rd-order solution ${\mathsf{\Theta }}_b^{\left[n=3\right]}$ to the rotation angle at the free end by IHAM vs. the point force curve is compared to the exact solutions ${\mathsf{\Theta }}_b$ in Figure 3. It can be stated that the improved homotopy analysis method has the optimal control parameters $\mathrm{ћ}$ for the convergence rate for an arbitrary $n$th-order approximation solution. For example, when $n=3$, the optimal control parameter is identified to be $\mathrm{ћ}=-0.5$.

4.4. Comparison between IHAM and HAM

The 30th-order approximation of ${\theta }_b$ with $n=3$, $\mathrm{ћ}=-0.5$ by HAM [44] can be calculated as:

$$\begin{array}{c}{\theta }_b=0.5\alpha -4.3435\times {10}^{-2}{\alpha }^3+6.8359\times {10}^{-3}{\alpha }^5-1.0902\times {10}^{-3}{\alpha }^7\hfill \\ +1.4940\times {10}^{-4}{\alpha }^9-1.6286\times {10}^{-5}{\alpha }^{11}+1.3403\times {10}^{-6}{\alpha }^{13}\hfill \\ -7.9070\times {10}^{-8}{\alpha }^{15}+3.1288\times {10}^{-9}{\alpha }^{17}-7.5167\times {10}^{-11}{\alpha }^{19}\hfill \\ +9.2764\times {10}^{-13}{\alpha }^{21}-4.2676\times {10}^{-15}{\alpha }^{23}+3.4351\times {10}^{-18}{\alpha }^{25}\hfill \\ -4.9207\times {10}^{-23}{\alpha }^{27}+1.4045\times {10}^{-32}{\alpha }^{29}\hfill \end{array}$$

(90)

When the relative error $\mathsf{\Delta }\le 1\%$, the convergent regions of the 30th-order approximation by HAM [44] are:

$$\alpha \in \left[0,4.3\right],\hspace{1em}{\mathsf{\Theta }}_b\in \left[0,0.74\right],$$

(91)

The 30th-order approximation of ${\theta }_b$ with $\epsilon =8/{\pi }^2$, $\mathrm{ћ}=-2.25$ by IHAM is

$$\begin{array}{c}{\theta }_b=0.5\alpha -4.1947\times {10}^{-2}{\alpha }^3+6.3195\times {10}^{-3}{\alpha }^5-1.0049\times {10}^{-3}\hfill \\ +1.4774\times {10}^{-4}{\alpha }^9-1.9106\times {10}^{-5}{\alpha }^{11}+2.1234\times {10}^{-6}{\alpha }^3\hfill \\ -2.0002\times {10}^{-7}{\alpha }^5+1.5807\times {10}^{-8}{\alpha }^{77}-1.0388\times {10}^{-9}{\alpha }^9\hfill \\ +5.619\times {10}^{-11}{\alpha }^{21}-2.5035\times {10}^{-12}{\alpha }^{23}+9.0941\times {10}^{-14}{\alpha }^3\hfill \\ -2.7016\times {10}^{-15}{\alpha }^{27}+6.5971\times {10}^{-17}{\alpha }^2-1.3356\times {10}^{-18}{\alpha }^{31}\hfill \\ +2.2653\times {10}^{-20}{\alpha }^{33}-3.2538\times {10}^{-22}{\alpha }^{35}+3.9973\times {10}^{-24}{\alpha }^{37}\hfill \\ -4.2352\times {10}^{-26}{\alpha }^{39}+3.8950\times {10}^{-8}{\alpha }^{41}-3.1220\times {10}^{-30}{\alpha }^{43}\hfill \\ +2.1842\times {10}^{-22}{\alpha }^{45}-1.3318\times {10}^{-34}{\alpha }^{47}+7.0399\times {10}^{-37}{\alpha }^{49}\hfill \\ -3.1913\times {10}^{-39}{\alpha }^{51}+1.2180\times {10}^{-41}{\alpha }^3-3.7894\times {10}^{-4}{\alpha }^{55}\hfill \\ +9.0212\times {10}^{-47}{\alpha }^5-1.4440\times {10}^{-98}{\alpha }^{59}+1.3901\times {10}^{-53}{\alpha }^{61}\hfill \end{array}$$

(92)

When the relative error $\mathsf{\Delta }\le 1\%$, the convergent regions of the 30th-order approximation by IHAM are:

$$\alpha \in \left[0,1.3\right],{\hspace{1em}\mathsf{\Theta }}_b\in \left[0,0.36\right],$$

(93)

Table 4. Thirtieth-order solutions ${\mathsf{\Theta }}_b^{\left[n=30\right]}$ by IHAM and HAM.

$\mathit{\alpha }$	Exact Solution	$\mathit{\epsilon }=0$ $\mathbf{ћ}=-1/10$ (HAM [48])		$\mathit{\epsilon }=8/{\mathit{\pi }}^2$ $\mathbf{ћ}=-2/25$ (IHAM)
	${\mathsf{\Theta }}_{\mathit{b}}$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$
0	0	0	0	0	0
0.3	0.094719	0.094757	0.0403%	0.094782	0.0664%
0.6	0.18509	0.18533	0.131%	0.18551	0.229%
0.9	0.26801	0.34295	0.217%	0.26911	0.411%
1.2	0.34206	0.29440	0.262%	0.29507	0.554%
1.5	0.40714	0.40824	0.268%	0.40974	0.637%
1.8	0.46394	0.46512	0.253%	0.46705	0.670%
2.1	0.51342	0.51462	0.233%	0.51685	0.668%
2.4	0.55661	0.55780	0.214%	0.56021	0.647%
2.7	0.59444	0.59563	0.200%	0.59809	0.615%
3	0.62772	0.62892	0.191%	0.63134	0.577%
3.3	0.65714	0.65838	0.189%	0.66067	0.538%
3.6	0.68327	0.68494	0.244%	0.68669	0.500%
3.9	0.70659	0.71070	0.584%	0.70987	0.464%
4.2	0.72749	0.73465	0.985%	0.73053	0.418%
4.5	0.74630	0.70350	5.73%	0.74859	0.307%
4.8	0.76329	0.28076	63.2%	0.76402	0.0959%
5.1	0.77870	—	—	0.78403	0.684%
5.4	0.79272	—	—	0.84973	7.19%
5.7	0.80552	—	—	—	—
6	0.81723	—	—	—	—
$\infty$	1	1	0	1	0

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 $\alpha \in \left[0,1.3\right]$ by the homotopy analysis method to $\mathrm{ћ}=-0.5$.

4.5. Homotopy–Páde Approximation

With application of the homotopy–Páde acceleration technique [26], the convergence region can be substantially enlarged. For example, when $\mathrm{ћ}=-30$, $\epsilon =0$, $\mathrm{ћ}=-1/10$, the $Pade\left[11,11\right]$ homotopy–Páde approximation of ${\theta }_b$ in Equation (92) by HAM [44] can be written as:

$${\theta }_b=\frac{\alpha }{2}\frac{f\left(\alpha \right)}{g\left(\alpha \right)}$$

(94a)

where

$$\begin{array}{c}f\left(\alpha \right)=1+0.19449{\alpha }^2+3.2631\times {10}^{-2}{\alpha }^4+2.2845\times {10}^{-3}{\alpha }^6\\ +1.8114\times {10}^{-4}{\alpha }^8+1.3556\times {10}^{-6}{\alpha }^{10}\end{array}$$

(94b)

$$\begin{array}{c}g\left(\alpha \right)=1+0.28136{\alpha }^2+4.3400\times {10}^{-2}{\alpha }^4+4.3884\times {10}^{-3}{\alpha }^6\\ +2.8367\times {10}^4{\alpha }^8+9.1351\times {10}^{-6}{\alpha }^{10}\end{array}$$

(94c)

When the relative error $\mathsf{\Delta }\le 1\%$, the convergent regions of the $Pade\left[11,11\right]$ homotopy–Páde approximation by HAM [44] are:

$$\alpha \in \left[0,4.8\right],{\hspace{1em}\mathsf{\Theta }}_b\in \left[0,0.76\right],$$

(95)

When $n=30,\epsilon =8/{\pi }^2,\mathrm{ћ}=-2/25$, the $Pade\left[15,15\right]$ homotopy–Páde approximation of ${\theta }_b$ in Equation (92) by IHAM is:

$${\theta }_b=\frac{\alpha }{2}\frac{f\left(\alpha \right)}{g\left(\alpha \right)}$$

(96a)

where

$$\begin{array}{c}f\left(\alpha \right)=1\&+0.27701{\alpha }^2+4.9969\times {10}^{-2}{\alpha }^4+5.2892\times {10}^{-3}{\alpha }^6+4.6124\times {10}^{-4}{\alpha }^8+2.2630\times {10}^{-5}{\alpha }^{10}\hfill \\ +8.9473\times {10}^{-7}{\alpha }^{12}+6.2843\times {10}^{-9}{\alpha }^{14}\hfill \end{array}$$

(96b)

$$\begin{array}{c}g\left(\alpha \right)=1+\&0.36091{\alpha }^2+6.7608\times {10}^{-2}{\alpha }^4+8.4095\times {10}^{-3}{\alpha }^6+7.4211\times {10}^{-4}{\alpha }^8+4.6051\times {10}^{-5}{\alpha }^{10}\hfill \\ +1.8474\times {10}^{-6}{\alpha }^{12}+3.6726\times {10}^{-8}{\alpha }^{14}\hfill \end{array}$$

(96c)

When the relative error $\mathsf{\Delta }\le 1\%$, the convergent region of the $Pade\left[15,15\right]$ homotopy–Páde approximation by IHAM are:

$$\alpha \in \left[0,7.6\right],{\hspace{1em}\mathsf{\Theta }}_b\in \left[0,0.87\right],$$

(97)

Table 5. Thirtieth homotopy–Páde approximation solutions ${\mathsf{\Theta }}_b^{\left[P\right]}$ by IHAM and HAM.

$\mathit{\alpha }$	Exact Solution	$\mathit{\epsilon }=0$ (HAM [48]) $\mathit{P}\mathit{a}\mathit{d}\mathit{e}\left[11,11\right]$		$\mathit{\epsilon }=8/{\mathit{\pi }}^2$ (IHAM) $\mathit{P}\mathit{a}\mathit{d}\mathit{e}\left[15,15\right]$
	${\mathsf{\Theta }}_{\mathit{b}}$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$	${\mathsf{\Theta }}_{\mathit{b}}$	$\mathsf{\Delta }$
0	0	0	0	0	0
0.5	0.15567	0.15583	0.0991%	0.15594	0.173%
1	0.29371	0.29440	0.237%	0.29507	0.463%
1.5	0.40714	0.40824	0.268%	0.40974	0.639%
2	0.49768	0.49887	0.239%	0.50102	0.671%
2.5	0.56977	0.57095	0.206%	0.57340	0.637%
3	0.62772	0.62875	0.165%	0.63134	0.577%
3.5	0.67489	0.67527	0.0565%	0.67834	0.511%
4	0.71380	0.71244	0.192%	0.71700	0.448%
4.5	0.74630	0.74176	0.608%	0.74920	0.389%
5	0.77373	0.76472	1.16%	0.77639	0.344%
5.5	0.79711	0.78279	1.80%	0.79971	0.326%
6	0.81723	0.79728	2.44%	0.82010	0.351%
6.5	0.83466	0.80932	3.04%	0.83834	0.441%
7	0.84986	0.81979	3.53%	0.85507	0.613%
7.5	0.86321	0.82935	3.92%	0.87080	0.879%
8	0.87499	0.83850	4.17%	0.88594	1.25%
8.5	0.88544	0.84758	4.27%	0.90079	1.73%
9	0.89475	0.85684	4.24%	0.91557	2.33%
9.5	0.90307	0.86644	4.06%	0.93046	3.03%
10	0.91055	0.87648	3.74%	0.94558	3.85%
$\infty$	1	1	0	1	0

compares the 30th homotopy–Páde approximation solutions by IHAM and HAM. From Table 5, the improved homotopy analysis method extends the convergence region from $\alpha \in \left[0,4.8\right],$ to ${\mathsf{\Theta }}_b\in \left[0,0.76\right]$ for the homotopy–Páde approximate solution against the original homotopy analysis method [44], which is a 58% extension.

Figure 4 compares the 30th-order iteration solution with the exact solution, from which the followings can be observed:

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

5. Discussion about the Relative Errors

In this section, the three iterations of $\alpha =2.0$ are taken as an example to discuss the effect of the control parameters $\epsilon$ and $\mathrm{ћ}$ on the relative errors.

After three iterations, the improved homotopy analysis method is adopted to solve Equations (68) and (69) whose solution in the following form, is the approximate of the original Equations (47) and (48)

$$\begin{array}{c}{\theta }_b=0.5\alpha +0.045833\mathrm{ћ}{\alpha }^3+0.091667\epsilon \mathrm{ћ}{\alpha }^3+0.13750\epsilon {\mathrm{ћ}}^2{\alpha }^3\hfill \\ +0.045833\epsilon {\mathrm{ћ}}^3{\alpha }^3+0.024400{\epsilon }^2{\mathrm{ћ}}^2{\alpha }^5+0.016267{\epsilon }^2{\mathrm{ћ}}^3{\alpha }^5+0.0018000{\epsilon }^3{\mathrm{ћ}}^3{\alpha }^7\hfill \end{array}$$

(98)

When $\alpha =2.0$, Equation (98) becomes:

$$\begin{gathered} {\theta }_b=1+0.36667\mathrm{ћ}+0.73333\epsilon \mathrm{ћ}+1.1000\epsilon {\mathrm{ћ}}^2+0.36667\epsilon {\mathrm{ћ}}^3+0.78079{\epsilon }^2{\mathrm{ћ}}^2 \\ +0.52053{\epsilon }^2{\mathrm{ћ}}^3+0.23040{\epsilon }^3{\mathrm{ћ}}^3 \end{gathered}$$

(99)

When $\alpha =2.0$, the exact solution to the ration angle at the free end is ${\mathsf{\Theta }}_b=0.49768$. When $\mathrm{ћ}=-0.45$, Equation (99) becomes:

$${\theta }_b=0.83500-0.14066\epsilon +0.11068{\epsilon }^2-0.020995{\epsilon }^3$$

(100)

According to Equation (100), the relative error curve with variable control parameter $\epsilon$ is a cubic curve. When the relative error is less than one percent, the convergence interval is $R_{\epsilon }=\left[0.5,1\right]$, which is the effective region of $\epsilon$. Figure 5 depicts the relative errors curve ${\mathsf{\Delta }}_{\mathsf{\Theta }}-\epsilon$ of original Equation (47) when, $\alpha =2.0,\mathrm{ћ}=-0.45$.

From Figure 5, the relative error is the minimum when $\epsilon =8/{\pi }^2$.

When $\epsilon =8/{\pi }^2$, Equation (99) becomes:

$${\theta }_b=1+0.96108\mathrm{ћ}+1.4046{\mathrm{ћ}}^2+0.76191{\mathrm{ћ}}^3$$

(101)

According to Equation (101), the relative error curve with variable control parameter $\mathrm{ћ}$ is also cubic. When the relative error is less than one percent, the convergence interval is $R_t=\left[-0.52,-0.4\right]$, which is the effective region of $\mathrm{ћ}$. Figure 6 shows the relative errors curve ${\mathsf{\Delta }}_{\mathsf{\Theta }}-\epsilon$ of original Equation (47) when $\alpha =20,\epsilon =8/{\pi }^2$.

From Figure 6, the relative error is the minimum when $\mathrm{ћ}=-0.45$.

6. 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 SOLUTION TO LARGE DEFORMATION OF CANTILEVER BEAM BY IMPROVED HOMOTOPY ANALYSIS METHOD II: VERTICAL AND HORIZONTAL DISPLACEMENTS”.