Abstract
In this paper we propose an accelerated convergence method, which is combined with the homotopy analysis method (HAM), to solve nonlinear problems. The HAM is applied to obtain approximate expressions. According to the numbers of terms in the approximations, some ratiocontrol parameters are introduced in the solution expressions. By solving simultaneous algebraic equations, all artificial parameters can be optimally identified, including the socalled convergencecontrol parameter $\hslash $. Twoexamples are given to illustrate the validity of the new method. Comparison with LP perturbation method and RungeKutta method reveals that the improved HAM is better than the standard HAM and applies especially to the problems with complicated nonlinear terms.
1. Introduction
There are various types of nonlinear problems in the engineering field, most of which are strongly nonlinear and do not have an analytical solution. To date, many powerful methods have been developed to solve nonlinear differential equations, such as the KrylovBogoliubovMitropolsky (KBM) [13], the method of harmonic balance [46], Adomian’s decomposition method (ADM) [7, 8] and the small parameter method [911]. However, these methods cannot provide us with a simple way to adjust and control the convergence region and rate of approximate solution series. The homotopy analysis method (HAM) [1215] was proposed in 1992 by Liao [16]. Based on the homotopy of topology, the advantage of the HAM is that it is independent of a small or large physical parameter and provides a way to control and ensure the convergence of approximate series. It has been proven that the HAM logically contains nonperturbations such as the Lyapunov artificial small parameter method, the ADM, the $\delta $expansion method, the variational iteration method and so on [17, 18]. Therefore, the HAM is more general.
However, to obtain a rather accurate homotopy approximation, the highorder deformation equation must be solved, which is usually complicated to compute. And in the standard HAM, the so called $\mathrm{\hslash}$curves must be plotted to find the optimal values of $\mathrm{\hslash}$ with which the approximations converge fastest. In this paper, we introduce an accelerated convergence method, by which the values of $\mathrm{\hslash}$ can be optimally recognized and the improved approximations converge faster than the standard homotopy approximations.
2. Basics of the HAM
One can transform a nonlinear problem into an infinite number of linear subproblems and obtain series solutions via the homotopy analysis method, whether the nonlinear problem contains small parameters or not. Consider the following general nonlinear system:
where $L$ is a linear operator, $N$ is a nonlinear operator, $g\left(r\right)$ is a known function, $u\left(r\right)$ is an unknown function to be determined and $r$ denotes temporal or spatial independent variables. Let ${u}_{0}\left(r\right)$ denote an initial guess of the exact solution for $u\left(r\right)$, $\mathrm{\hslash}\ne 0$ denote an auxiliary convergencecontrol parameter, and $H\left(r\right)\ne 0$ denote an auxiliary function. Then, the homotopy function is constructed as:
where $q\in \left[\mathrm{0,1}\right]$ is an embedding parameter, $\mathcal{L}$ is an auxiliary linear operator with the property $\mathcal{L}\left(0\right)=0$, and $\mathrm{\aleph}$ is a nonlinear operator defined as:
Then, enforcing:
the zeroorder deformation equation would be constructed as:
When $q=\text{0}$ and $q=\text{1}$, we have respectively:
when $q=\text{0}$, and:
when $q=\text{1}$.
It is easy to show that from Eq. (6) and Eq. (7):
Thus, as the embedding parameter $q$ increases from 0 to 1, $\mathrm{\Phi}(r;q)$ varies continuously from the initial guess ${u}_{0}\left(r\right)$ to the exact solution $u\left(r\right)$ of the original Eq. (1). Differentiating Eq. (5) $n$ times with respect to the embedding parameter $q$ and then setting $q=\text{0}$ and finally dividing them by $n!$, one can obtain the socalled $n$thorder deformation equation:
where:
for simplicity, and:
Because $\mathrm{\Phi}(r;q)$ depends on the embedding parameter $q\in \left[\mathrm{0,1}\right]$, it can be expanded into a power series of $q$ by using the Tylor theorem as:
where:
${u}_{m}\left(r\right)$ can be found by solving Eq. (10), which is a linear equation. If the initial guess ${u}_{0}\left(r\right)$, auxiliary linear operator $\mathcal{L}$, convergencecontrol parameter $\mathrm{\hslash}$ and auxiliary function $H\left(r\right)$ are properly chosen, the series Eq. (14) is convergent to $\mathrm{\Phi}(r;q)$ at $q=\text{1}$, and one can have, according to Eq. (9) and Eq. (14), the homotopyseries solution:
and the $n$thorder homotopyapproximation:
3. The idea of the accelerated convergence method
Only the series solutions can be obtained for most nonlinear problems. It makes sense to ensure and accelerate the convergence of the series. Once the base functions that express the solution are chosen, the proportion of the terms of the solution is pivotal to the convergence rate. We try to introduce some ratiocontrol parameters in the approximate solution expressions, combined with the original governing equations and the given conditions, to control the proportion between the terms. If there are $n+1$ nonconstant terms in the approximate expressions, $n$ ratiocontrol parameters should be introduced as multipliers of the first $n$ terms. Assume that the $x$thorder homotopy approximation of Eq. (1) can be expressed as:
Then we introduce $n$ ratiocontrol parameters in Eq. (18) as follows:
For initial value or boundary value problems, the values of $u\left(s\right)$, $u\text{'}\left(s\right)$, …, and ${u}^{\left(n\right)}\left(s\right)$ can be obtained by continuous derivation of Eq. (1) if the given conditions is at $r=s$, then an equation can be written as follows:
Note that the numbers of equations should be equal to the numbers of unknowns in Eq. (20) including the convergencecontrol parameter $\mathrm{\hslash}$. Sometimes Eq. (20) consists of nonlinear algebraic equations, thus there may be multiple solutions. Therefore, to choose the proper values of parameters, we can use the squared residual:
where, $v$ denotes the number of the solutions of Eq. (20) and:
Then, one can obtain the optimal values of parameters for which ${E}_{x}({a}_{1,v},{a}_{2,v},\cdots ,{a}_{n,v},{\mathrm{\hslash}}_{v})$ is the minimum.
Two examples are evaluated in the following to illustrate the specific method.
4. Illustrative examples
4.1. Duffing equation
First, consider the duffing equation:
where the dot denotes the derivative with respect to $t$. The base functions and the initial approximation of the solution are chosen respectively as:
where $\omega $ is the frequency to be determined and:
The auxiliary linear operator can be chosen as:
and the nonzero auxiliary function $H\left(t\right)$ as 1 for simplicity.
By using the homotopy analysis method, one can easily obtain the secondorder homotopy approximations of Eq. (23) as:
where:
There are three nonconstant terms in the secondorder approximation ${u}_{2}\left(t\right)$. Thus, we introduce two socalled ratiocontrol parameters $\mu $ and $\lambda $ into Eq. (28) as:
One can easily obtain the expressions of ${\ddot{\stackrel{~}{u}}}_{2}\left(t\right)$, ${\stackrel{~}{u}}_{2}^{\left(4\right)}\left(t\right)$ and ${\stackrel{~}{u}}_{2}^{\left(6\right)}\left(t\right)$ by differentiating Eq. (29) continuously with respect to $t$, respectively and the exact values of $u\left(0\right)$, $\ddot{u}\left(0\right)$, ${u}^{\left(4\right)}\left(0\right)$ and ${u}^{\left(6\right)}\left(0\right)$ from Eq. (23). Then, simultaneous equations can be written as follows:
In order to choose the proper values of unknowns, for the sake of computational efficiency, we calculate the discrete squared residual:
where $m$ denotes the order of approximations, $N$ is an integer and:
Obviously, Eq. (31) is a good approximation of Eq. (21) for large enough $N$. In this example, we use $N=$ 100. The values of the unknows $\beta $, $\mu $, $\lambda $ and $\mathrm{\hslash}$ can be obtained by solving Eq. (30) and using Eq. (31). Then the improved frequency ${\stackrel{~}{\omega}}_{2}$ can be obtained from Eq. (27). The second order standard homotopy approximations can be easily obtained by plotting the so called $\mathrm{\hslash}$curves for different cases. Table 1 gives the values of parameter $\mathrm{\hslash}$ which minimizes the discrete squared residual in the standard HAM (SHAM). To illustrate the validity of the improved homotopy analysis method (IHAM) for this example, comparison of frequency is given in Table 2.
Table 1The minimum of discrete squared residual in different cases for example 1
$\epsilon =1$  
A  0.1  1.0  10.0  100.0  1000.0 
$\hslash $  –0.9826  –0.4938  –0.01015  –1.026×10^{4}  –1.027×10^{5} 
SHAM  4.824×10^{12}  1.302×10^{2}  5.723×104  5.851×10^{10}  5.852×10^{16} 
IHAM  1.383×10^{16}  2.054×10^{3}  7.476×10^{4}  7.898×10^{10}  7.902×10^{16} 
Table 2Comparison of ω corresponding to various parameters of system for example 1
$\epsilon =1$  
A  0.1  1.0  10.0  100.0  1000.0 
${\omega}_{2}$  1.00374337  1.32387414  8.74071293  86.84060261  868.34896697 
${\stackrel{~}{\omega}}_{2}$  1.00374184  1.31893874  8.65607864  85.99897789  859.93345617 
Exact values  1.00374184  1.31777607  8.53358619  84.72747994  847.21370197 
4.2. Simple pendulum attached to a rotating rigid frame[12]
The governing equation is:
subject to the initial conditions:
To apply the HAM to solve the above problem, the following approximation is used:
as a result, we obtain the following approximate differential equation:
Obviously, the more terms of the Maclaurin series considered, the more efficient the solution series obtained but the more complicated the problem will be. To illustrate the validity and simplicity of the accelerated convergence method, we only use the first two terms of the Maclaurin series to approximate the original Eq. (34) as a duffinglike Eq. (37). The base functions and the initial approximation of the solution are chosen as Eq. (24) and Eq. (25), respectively. The second order approximate expressions of Eq. (37) can be easily obtained by the HAM. There are three nonconstant terms in the secondorder approximate expression ${\theta}_{2}\left(t\right)$. Thus, the improved expression ${\stackrel{~}{\theta}}_{2}\left(t\right)$ should contain two ratiocontrol parameters like Eq. (29).
The second order homotopy approximate frequency of Eq. (37) reads:
Note that, to establish the simultaneous equations like Eq. (30) in this example, one should obtain the exact values of $\theta \left(0\right)$, $\ddot{\theta}\left(0\right)$, ${\theta}^{\left(4\right)}\left(0\right)$ and ${\theta}^{\left(6\right)}\left(0\right)$ by differentiating the original governing Eq. (34) continuously with respect to $t$ and using the initial conditions.
We also obtain the values of convergencecontrol parameter $\mathrm{\hslash}$ which minimizes the discrete squared residual in the SHAM by plotting the $\mathrm{\hslash}$curves and compare the minimum of discrete squared residual in the SHAM and the IHAM, as shown in Table 3. To illustrate the validity of the accelerated convergence method for this example, comparison of frequency is given in Table 4.
Table 3The minimum of discrete squared residual in different cases for example 2
$\omega =$ 1, $\epsilon =$0.5, $\varphi =\text{0.5}\pi $  $\omega =$1,$\epsilon =$0.5, $\varphi =\text{0.7}\pi $  $\omega =$1,$\epsilon =$0.9, $\varphi =\text{0.5}\pi $  $\omega =$1,$\epsilon =$0.9, $\varphi =\text{0.7}\pi $  $\omega =$2,$\epsilon =$0.9, $\varphi =\text{0.7}\pi $  
$\hslash $  –0.503  –0.189  –0.219  –0.090  –0.032 
SHAM  128.28  367.91  138.00  756.54  12048.3 
IHAM  0.014  0.722  0.003  2.780  44.475 
Table 4Comparison of α corresponding to various parameters of system for example 2
$\omega =$ 1, $\epsilon =$0.5, $\varphi =0.5\pi $  $\omega =$1,$\epsilon =$0.5, $\varphi =0.7\pi $  $\omega =$1,$\epsilon =$0.9, $\varphi =0.5\pi $  $\omega =$1,$\epsilon =$0.9, $\varphi =0.7\pi $  $\omega =$2,$\epsilon =$0.9, $\varphi =0.7\pi $  
${\alpha}_{2}$  0.90006070  1.05130991  0.94966034  1.29215650  2.58072027 
${\stackrel{~}{\alpha}}_{2}$  0.79728921  0.75460788  0.74192670  0.79527954  1.59055907 
exact values  0.79432027  0.74093414  0.74305823  0.76644529  1.53289057 
Fig. 1Comparison of the exact solution with the analytic approximations in the case of ω= 1, ε= 0.5, ϕ=0.5π for example 2. Symbols: numerical solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting ℏcurves
Fig. 2Comparison of the exact solution with the analytic approximations in the case of ω= 2, ε= 0.9, ϕ=0.7π for example 2. Symbols: numerical solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting ℏcurves
5. Comparison with LP perturbation method
LindstedtPoincaré perturbation method [19, 20] (LP perturbation method) is a powerful analytical technique to weak nonlinear problems. Eq. (34) have a periodic solution. In the idea of LP perturbation, one should express the solution $u\left(t\right)$ and the frequency squared ${\omega}^{2}$ in power series with respect to $\epsilon $:
where, $\left\epsilon \right\ll \text{1}$. By substituting Eq. (39) and Eq. (40) into Eq. (23), one obtain:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\epsilon ({u}_{0}+\epsilon {u}_{1}+{\epsilon}^{2}{u}_{2}+\cdots {)}^{3}=0,$
By equating coefficients of like powers of $\epsilon $, one can obtain the set of recursive linear differential equations:
and so on. By solving above equations and eliminating the secular terms, we can obtain three terms order LP approximation as:
Fig. 3 and Fig. 4 give the solution comparisons of LP perturbation method with SHAM and IHAM in different weak nonlinear cases.
Fig. 3Comparison of the analytic approximations in the case of ε= 0.1, A= 5 for example 1. Symbols: LP perturbation solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting ℏcurves
Fig. 4Comparison of the analytic approximations in the case of ε= 0.01, A= 10 for example 1. Symbols: LP perturbation solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting ℏcurves
6. Results and discussion
The standard homotopy analysis method (SHAM) gives us a way to control the ratio and region of approximations for nonlinear problems by choosing the values of $\mathrm{\hslash}$, thus in general it can give better solutions than other analytical methods. Without plotting $\mathrm{\hslash}$curves, the same and even more accurate results can be obtained by improved homotopy analysis method (IHAM) in various cases, as shown in Table 1. Table 2 shows that the accelerated convergence method can give more accurate values of frequencies than $\mathrm{\hslash}$curves. We also applied the IHAM to a problem with complicated nonlinear terms and obtained better results than ones by plotting the $\mathrm{\hslash}$curves in different cases, as shown in Table 4, Table 5, Fig. 1 and Fig. 2. This successful application means that the problems with complicated nonlinear terms can be reduced to simple ones. One can only analyze simplified duffinglike equations by the HAM and obtain greatly accurate results of the original problems by the accelerated convergence method. For weak nonlinear problems, LP perturbation can give accurate enough approximate solution. Comparison with LP perturbation reveals that IHAM is of high accuracy for various cases, as shown in Fig. 3 and Fig. 4.
7. Conclusions
In this paper, we proposed an accelerated convergence method and improved the HAM. Faster convergent results can be obtained by our method than plotting the $\mathrm{\hslash}$curves. The idea of the method is to control the proportion between the homotopy approximation terms by introducing some ratiocontrol parameters in the approximate solution expressions. All artificial parameters including the so called convergencecontrol parameter $\mathrm{\hslash}$ can be optimally identified by solving algebraic equations. The current work has elucidated the validity and justifications of this combination.
References

Krylov N. M., Bogoli︠u︡bov N. M. Introduction to NonLinear Mechanics. Princeton University Press, 1943.

Kakutani T., Sugimoto N. KrylovBogoliubovMitropolsky method for nonlinear wave modulation. Physics of Fluids, Vol. 17, Issue 8, 1974, p. 16171625.

Yamgoué S. B., Kofané T. C. Application of the KrylovBogoliubovMitropolsky method to weakly damped strongly nonlinear planar Hamiltonian systems. International Journal of NonLinear Mechanics, Vol. 42, Issue 10, 2007, p. 12401247.

Hu H. Solution of a quadratic nonlinear oscillator by the method of harmonic balance. Journal of Sound and Vibration, Vol. 293, Issue 1, 2006, p. 462468.

Ghadimi M., Kaliji H. Application of the harmonic balance method on nonlinear equations. World Applied Sciences Journal, Vol. 22, Issue 4, 2013, p. 532537.

Ju P., Xue X. Global residue harmonic balance method for largeamplitude oscillations of a nonlinear system. Applied Mathematical Modelling, Vol. 39, Issue 2, 2015, p. 449454.

ElSayed A., Behiry S., Raslan W. Adomian’s decomposition method for solving an intermediate fractional advectiondispersion equation. Computers and Mathematics with Applications, Vol. 59, Issue 5, 2010, p. 17591765.

Ebaid A., Aljoufi M. D., Wazwaz A. An advanced study on the solution of nanofluid flow problems via Adomian’s method. Applied Mathematics Letters, Vol. 46, 2015, p. 117122.

Andrianov I., Awrejcewicz J., Ivankov A. Artificial small parameter methodsolving mixed boundary value problems. Mathematical Problems in Engineering, Vol. 3, Issue 2005, 2005, p. 325340.

Andrianov I., Danishevs’Kyy V., Awrejcewicz J. An artificial small perturbation parameter and nonlinear plate vibrations. Journal of Sound and Vibration, Vol. 283, Issue 3, 2005, p. 561571.

Ramos J. An artificial parameterdecomposition method for nonlinear oscillators: Applications to oscillators with odd nonlinearities. Journal of Sound and Vibration, Vol. 307, Issue 1, 2007, p. 312329.

Liao S. J., Chwang A. Application of homotopy analysis method in nonlinear oscillations. Journal of Applied Mechanics, Vol. 65, Issue 4, 1998, p. 914922.

Liao S. J., Cheung K. F. Homotopy analysis of nonlinear progressive waves in deep water. Journal of Engineering Mathematics, Vol. 45, Issue 2, 2003, p. 105116.

Liao S. J. On the homotopy analysis method for nonlinear problems. Applied Mathematics and Computation, Vol. 147, Issue 2, 2004, p. 499513.

Park S. H., Kim J. H. Homotopy analysis method for option pricing under stochastic volatility. Applied Mathematics Letters, Vol. 24, Issue 10, 2011, p. 17401744.

Liao S. J. The Proposed Homotopy Analysis Technique for the Solution of Nonlinear Problems. Ph.D. Thesis, Shanghai Jiao Tong University, 1992.

Liao S. J. Beyond Perturbation: Introduction to the Homotopy Analysis Method. CRC Press, 2003.

Van Gorder R. A. The variational iteration method is a special case of the homotopy analysis method. Applied Mathematics Letters, Vol. 45, 2015, p. 8185.

Nayfeh A. H. Perturbation Methods. John Wiley and Sons, 2008.

Marinca V., Herisanu N. Nonlinear Dynamical Systems in Engineering. Springer Berlin Heidelberg, 2012, p. 929.
About this article
The authors appreciate the financial support of the National Natural Science Foundation of China (Grant Number 51205167) and that of the Natural Science Foundation of Jiangsu Province, China (Grant Number BK20151128), the financial support from State Key Laboratory for Strength and Vibration of Mechanical Structures (Grant Number SV2015KF11).