An improved homotopy analysis method with accelerated convergence for nonlinear problems

Hong-wei Li1 , Jun Wang2 , Li-xin Lu3 , Zhi-wei Wang4

1, 2, 3Department of Packaging Engineering, Jiangnan University, Wuxi 214122, China

2State Key Laboratory for Strength and Vibration of Mechanical Structures, Xian 710049, China

4Packaging Engineering Institute, Jinan University, Zhuhai 519070, China

2Corresponding author

Journal of Vibroengineering, Vol. 18, Issue 7, 2016, p. 4756-4765. https://doi.org/10.21595/jve.2016.16808
Received 8 January 2016; received in revised form 14 July 2016; accepted 18 August 2016; published 15 November 2016

Copyright © 2016 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Table of Contents Acknowledgements References
Cite this article
Views 131
Reads 73
Downloads 1212
WoS Core Citations 0
CrossRef Citations 0
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 ratio-control parameters are introduced in the solution expressions. By solving simultaneous algebraic equations, all artificial parameters can be optimally identified, including the so-called convergence-control parameter . Twoexamples are given to illustrate the validity of the new method. Comparison with L-P perturbation method and Runge-Kutta method reveals that the improved HAM is better than the standard HAM and applies especially to the problems with complicated nonlinear terms.

Keywords: nonlinear differential equations, homotopy analysis method, analytical approximation, convergence acceleration.

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 Krylov-Bogoliubov-Mitropolsky (KBM) [1-3], the method of harmonic balance [4-6], Adomian’s decomposition method (ADM) [7, 8] and the small parameter method [9-11]. 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) [12-15] 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 non-perturbations such as the Lyapunov artificial small parameter method, the ADM, the δ-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 high-order deformation equation must be solved, which is usually complicated to compute. And in the standard HAM, the so called -curves must be plotted to find the optimal values of with which the approximations converge fastest. In this paper, we introduce an accelerated convergence method, by which the values of 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 sub-problems and obtain series solutions via the homotopy analysis method, whether the nonlinear problem contains small parameters or not. Consider the following general nonlinear system:

(1)
L u r + N u r = g r ,

where L is a linear operator, N is a nonlinear operator, g(r) is a known function, u(r) is an unknown function to be determined and r denotes temporal or spatial independent variables. Let u0(r) denote an initial guess of the exact solution for u(r), 0 denote an auxiliary convergence-control parameter, and H(r)0 denote an auxiliary function. Then, the homotopy function is constructed as:

(2)
H Φ r ; q , u 0 r , H r , h , q = 1 - q L Φ r ; q - u 0 r - q H r Φ r ; q ,

where q[0,1] is an embedding parameter, L is an auxiliary linear operator with the property L(0)=0, and is a nonlinear operator defined as:

(3)
Φ r ; q = L Φ r ; q + N Φ r ; q - g r .

Then, enforcing:

(4)
H [ Φ ( r ; q ) , u 0 ( r ) , H ( r ) , h , q ] = 0 ,

the zero-order deformation equation would be constructed as:

(5)
1 - q L Φ r ; q - u 0 r = q H r Φ r ; q .

When q=0 and q=1, we have respectively:

(6)
L [ Φ ( r ; 0 ) - u 0 ( r ) ] = 0 ,

when q=0, and:

(7)
[ Φ ( r ; 1 ) ] = 0 ,

when q=1.

It is easy to show that from Eq. (6) and Eq. (7):

(8)
Φ r ; 0 = u 0 r ,
(9)
Φ r ; 1 = u r .

Thus, as the embedding parameter q increases from 0 to 1, Φ(r;q) varies continuously from the initial guess u0(r) to the exact solution u(r) of the original Eq. (1). Differentiating Eq. (5) n times with respect to the embedding parameter q and then setting q=0 and finally dividing them by n!, one can obtain the so-called nth-order deformation equation:

(10)
L u n r - χ n u n - 1 r = H r R n u n - 1 ,

where:

(11)
u m = u 0 ( r ) , u 1 ( r ) , u 2 ( r ) , , u m ( r ) ,

for simplicity, and:

(12)
R n u n - 1 = 1 n - 1 ! n - 1 Φ r ; q q n - 1   q = 0 ,
(13)
χ n = 0 ,       n 1 , 1 ,       otherwise .

Because Φ(r;q) depends on the embedding parameter q[0,1], it can be expanded into a power series of q by using the Tylor theorem as:

(14)
Φ r ; q ~ u 0 r + m = 1 + u m r q m ,

where:

(15)
u m r = 1 m ! m Φ r ; q q m   q = 0 .

u m ( r ) can be found by solving Eq. (10), which is a linear equation. If the initial guess u0(r), auxiliary linear operator L, convergence-control parameter and auxiliary function H(r) are properly chosen, the series Eq. (14) is convergent to Φ(r;q) at q=1, and one can have, according to Eq. (9) and Eq. (14), the homotopy-series solution:

(16)
u r = u 0 r + m = 1 + u m r ,

and the nth-order homotopy-approximation:

(17)
u r u 0 r + m = 1 n u m r .

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 ratio-control 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 non-constant terms in the approximate expressions, n ratio-control parameters should be introduced as multipliers of the first n terms. Assume that the xth-order homotopy approximation of Eq. (1) can be expressed as:

(18)
u x r u x , 1 r + u x , 2 r + + u x , n r + u x , n + 1 r + C .

Then we introduce n ratio-control parameters in Eq. (18) as follows:

(19)
u ~ x r a 1 u x , 1 r + a 2 u x , 2 r + + a n u x , n r + u x , n + 1 r + C .

For initial value or boundary value problems, the values of u(s), u'(s), …, and u(n)(s) 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:

(20)
a 1 u x , 1 s + a 2 u x , 2 s + + a n u x , n s + u x , n + 1 s + C = u s , a 1 u x , 1 ' s + a 2 u x , 2 ' s + + a n u x , n ' s + u x , n + 1 ' s = u ' s , a 1 u n x , 1 s + a 2 u n x , 2 s + + a n u n x , n s + u n x , n + 1 s = u n s .

Note that the numbers of equations should be equal to the numbers of unknowns in Eq. (20) including the convergence-control parameter . 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:

(21)
E x a 1 , v , a 2 , v , , a n , v , v = 0 r R x r ; a 1 , v , a 2 , v , , a n , v , v 2 d r ,

where, v denotes the number of the solutions of Eq. (20) and:

(22)
R x r ; a 1 , v , a 2 , v , , a n , v , v = L u ~ x r + N u ~ x r - g r .

Then, one can obtain the optimal values of parameters for which Ex(a1,v,a2,v,,an,v,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:

(23)
u ¨ + u + ε u 3 = 0 , u ( 0 ) = A ,         u ˙ ( 0 ) = 0 ,

where the dot denotes the derivative with respect to t. The base functions and the initial approximation of the solution are chosen respectively as:

(24)
c o s 2 m + 1 ω t m = 0 ,   1 ,   2 ,   ,

where ω is the frequency to be determined and:

(25)
u 0 t = β c o s ω t .

The auxiliary linear operator can be chosen as:

(26)
L u = u ¨ + u ,

and the non-zero auxiliary function H(t) as 1 for simplicity.

By using the homotopy analysis method, one can easily obtain the second-order homotopy approximations of Eq. (23) as:

(27)
ω 2 = 1 + 3 ε β 2 4 + ε 2 β 4 64 + 3 2 ε 2 β 4 128 + 39 2 ε 2 β 6 2048 ,
(28)
u 2 t = c 1 c o s ω 2 t + c 2 c o s 3 ω 2 t + c 3 c o s 5 ω 2 t ,

where:

c 1 = β + ε β 3 16 + 2 ε β 3 32 + 23 2 ε 2 β 5 1024 ,         c 2 = - ε β 3 16 + 2 ε β 3 32 + 3 2 ε 2 β 5 128 ,         c 3 = 2 ε 2 β 5 1024 .

There are three non-constant terms in the second-order approximation u2(t). Thus, we introduce two so-called ratio-control parameters μ and λ into Eq. (28) as:

(29)
u ~ 2 t = μ c 1 c o s ω 2 t + λ c 2 c o s 3 ω 2 t + c 3 c o s 5 ω 2 t .

One can easily obtain the expressions of u~¨2(t), u~2(4)(t) and u~2(6)(t) by differentiating Eq. (29) continuously with respect to t, respectively and the exact values of u(0), u¨(0), u(4)(0) and u(6)(0) from Eq. (23). Then, simultaneous equations can be written as follows:

(30)
μ c 1 + λ c 2 + c 3 = A , ω 2 2 μ c 1 + 9 λ c 2 + 25 c 3 = A + ε A 3 , ω 2 4 μ c 1 + 81 λ c 2 + 625 c 3 = A + 4 ε A 3 + 3 ε 2 A 5 , ω 2 6 μ c 1 + 729 λ c 2 + 15625 c 3 = A + 25 ε A 2 + 51 ε 2 A 4 + 27 ε 3 A 6 .

In order to choose the proper values of unknowns, for the sake of computational efficiency, we calculate the discrete squared residual:

(31)
E m ,   β ,   μ k = 0 N [ R m ( τ k ;   ,   β ,   μ ) ] 2 ,

where m denotes the order of approximations, N is an integer and:

(32)
τ k = 2 k π N ,
(33)
R m τ k ;   ,   β ,   μ = ω 2 u ~ ¨ τ k + u ~ τ k + ε u ~ 3 τ k .

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 β, μ, λ and can be obtained by solving Eq. (30) and using Eq. (31). Then the improved frequency ω~2 can be obtained from Eq. (27). The second order standard homotopy approximations can be easily obtained by plotting the so called -curves for different cases. Table 1 gives the values of parameter 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 1. The minimum of discrete squared residual in different cases for example 1

ε = 1
A
0.1
1.0
10.0
100.0
1000.0
–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×1010
5.852×1016
IHAM
1.383×10-16
2.054×10-3
7.476×104
7.898×1010
7.902×1016

Table 2. Comparison of ω corresponding to various parameters of system for example 1

ε = 1
A
0.1
1.0
10.0
100.0
1000.0
ω 2
1.00374337
1.32387414
8.74071293
86.84060261
868.34896697
ω ~ 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:

(34)
θ ¨ + ω 2 ( 1 - ε c o s θ ) s i n θ = 0 ,

subject to the initial conditions:

(35)
θ 0 = ϕ ,           θ ˙ ( 0 ) = 0 .

To apply the HAM to solve the above problem, the following approximation is used:

(36)
ω 2 1 - ε c o s θ s i n θ - ω 2 ε - 1 θ + ω 2 2 ε 3 - 1 6 θ 3 ,

as a result, we obtain the following approximate differential equation:

(37)
θ ¨ - ω 2 ( ε - 1 ) θ + ω 2 2 ε 3 - 1 6 θ 3 = 0 .

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 duffing-like 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 non-constant terms in the second-order approximate expression θ2(t). Thus, the improved expression θ~2(t) should contain two ratio-control parameters like Eq. (29).

The second order homotopy approximate frequency of Eq. (37) reads:

(38)
α 2 = ( 1 - ε ) ω 2 + 3 β 2 ω 2 4 2 ε 3 - 1 6 + β 4 ω 4 64 2 ε 3 - 1 6 2       + 3 2 β 4 ω 4 128 2 ε 3 - 1 6 2 + 39 2 β 6 ω 4 2048 2 ε 3 - 1 6 2 .

Note that, to establish the simultaneous equations like Eq. (30) in this example, one should obtain the exact values of θ(0), θ¨(0), θ(4)(0) and θ(6)(0) by differentiating the original governing Eq. (34) continuously with respect to t and using the initial conditions.

We also obtain the values of convergence-control parameter which minimizes the discrete squared residual in the SHAM by plotting the -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 3. The minimum of discrete squared residual in different cases for example 2

ω = 1, ε= 0.5,
ϕ = 0.5 π
ω =   1, ε= 0.5,
ϕ = 0.7 π
ω =   1, ε= 0.9,
ϕ = 0.5 π
ω =   1, ε= 0.9,
ϕ = 0.7 π
ω =   2, ε= 0.9,
ϕ = 0.7 π
–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 4. Comparison of α corresponding to various parameters of system for example 2

ω = 1, ε= 0.5,
ϕ = 0.5 π
ω =   1, ε= 0.5,
ϕ = 0.7 π
ω =   1, ε= 0.9,
ϕ = 0.5 π
ω =   1, ε= 0.9,
ϕ = 0.7 π
ω =   2, ε= 0.9,
ϕ = 0.7 π
α 2
0.90006070
1.05130991
0.94966034
1.29215650
2.58072027
α ~ 2
0.79728921
0.75460788
0.74192670
0.79527954
1.59055907
exact values
0.79432027
0.74093414
0.74305823
0.76644529
1.53289057

Fig. 1. Comparison 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

 Comparison 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. 2. Comparison 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

 Comparison 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 L-P perturbation method

Lindstedt-Poincaré perturbation method [19, 20] (L-P perturbation method) is a powerful analytical technique to weak nonlinear problems. Eq. (34) have a periodic solution. In the idea of L-P perturbation, one should express the solution u(t) and the frequency squared ω2 in power series with respect to ε:

(39)
u t = u 0 t + ε u 1 t + ε 2 u 2 t + ,
(40)
ω 2 = 1 + ε ω 1 + ε 2 ω 2 + ,

where, ε1. By substituting Eq. (39) and Eq. (40) into Eq. (23), one obtain:

(41)
u ¨ 0 + ε u ¨ 1 + ε 2 u ¨ 2 + ( ω 2 - ε ω 1 - ε ω 2 - ) ( u 0 + ε u 1 + ε 2 u 2 + )
            + ε ( u 0 + ε u 1 + ε 2 u 2 + ) 3 = 0 ,
(42)
u 0 0 + ε u 1 0 + ε 2 u 2 0 + = A , u ˙ 0 ( 0 ) + ε u ˙ 1 ( 0 ) + ε 2 u ˙ 2 ( 0 ) + = 0 .

By equating coefficients of like powers of ε, one can obtain the set of recursive linear differential equations:

(43)
ε 0 :   u ¨ 0 + ω 2 u 0 = 0 , u 0 0 = A ,             u ˙ 0 ( 0 ) = 0 ,
(44)
ε 1 :   u ¨ 1 + ω 2 u 1 - ω 1 u 0 + u 0 3 = 0 , u 1 0 = 0 ,           u ˙ 1 ( 0 ) = 0 ,
(45)
ε 2 :   u ¨ 2 + ω 2 u 2 - ω 1 u 1 - ω 2 u 0 + 3 u 0 2 u 1 = 0 , u 2 0 = 0 ,           u ˙ 2 ( 0 ) = 0 ,

and so on. By solving above equations and eliminating the secular terms, we can obtain three terms order L-P approximation as:

(46)
ω = 8 + 6 ε A 2 + 64 + 96 ε A 2 + 30 ε 2 A 4 4 ,
(47)
u = A - ε A 3 32 ω 2 - ε 2 A 5 1024 ω 4 c o s ω t + ε A 3 32 ω 2 c o s 3 ω t + ε 2 A 5 1024 ω 4 c o s 5 ω t .

Fig. 3 and Fig. 4 give the solution comparisons of L-P perturbation method with SHAM and IHAM in different weak nonlinear cases.

Fig. 3. Comparison of the analytic approximations in the case of ε= 0.1, A= 5 for example 1. Symbols: L-P perturbation solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting -curves

 Comparison of the analytic approximations in the case of ε= 0.1, A= 5 for example 1.  Symbols: L-P perturbation solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting ℏ-curves

Fig. 4. Comparison of the analytic approximations in the case of ε= 0.01, A= 10 for example 1. Symbols: L-P perturbation solution; solid line: improved homotopy approximation by accelerated convergence; dashed line: standard homotopy approximation by plotting -curves

 Comparison of the analytic approximations in the case of ε= 0.01, A= 10 for example 1.  Symbols: L-P 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 , thus in general it can give better solutions than other analytical methods. Without plotting -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 -curves. We also applied the IHAM to a problem with complicated nonlinear terms and obtained better results than ones by plotting the -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 duffing-like equations by the HAM and obtain greatly accurate results of the original problems by the accelerated convergence method. For weak nonlinear problems, L-P perturbation can give accurate enough approximate solution. Comparison with L-P 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 -curves. The idea of the method is to control the proportion between the homotopy approximation terms by introducing some ratio-control parameters in the approximate solution expressions. All artificial parameters including the so called convergence-control parameter can be optimally identified by solving algebraic equations. The current work has elucidated the validity and justifications of this combination.

Acknowledgements

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 SV2015-KF-11).

References

  1. Krylov N. M., Bogoli︠u︡bov N. M. Introduction to Non-Linear Mechanics. Princeton University Press, 1943. [Search CrossRef]
  2. Kakutani T., Sugimoto N. Krylov-Bogoliubov-Mitropolsky method for nonlinear wave modulation. Physics of Fluids, Vol. 17, Issue 8, 1974, p. 1617-1625. [Publisher]
  3. Yamgoué S. B., Kofané T. C. Application of the Krylov-Bogoliubov-Mitropolsky method to weakly damped strongly non-linear planar Hamiltonian systems. International Journal of Non-Linear Mechanics, Vol. 42, Issue 10, 2007, p. 1240-1247. [Publisher]
  4. 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. 462-468. [Publisher]
  5. Ghadimi M., Kaliji H. Application of the harmonic balance method on nonlinear equations. World Applied Sciences Journal, Vol. 22, Issue 4, 2013, p. 532-537. [Search CrossRef]
  6. Ju P., Xue X. Global residue harmonic balance method for large-amplitude oscillations of a nonlinear system. Applied Mathematical Modelling, Vol. 39, Issue 2, 2015, p. 449-454. [Publisher]
  7. El-Sayed A., Behiry S., Raslan W. Adomian’s decomposition method for solving an intermediate fractional advection-dispersion equation. Computers and Mathematics with Applications, Vol. 59, Issue 5, 2010, p. 1759-1765. [Publisher]
  8. 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. 117-122. [Publisher]
  9. Andrianov I., Awrejcewicz J., Ivankov A. Artificial small parameter method-solving mixed boundary value problems. Mathematical Problems in Engineering, Vol. 3, Issue 2005, 2005, p. 325-340. [Search CrossRef]
  10. 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. 561-571. [Publisher]
  11. Ramos J. An artificial parameter-decomposition method for nonlinear oscillators: Applications to oscillators with odd nonlinearities. Journal of Sound and Vibration, Vol. 307, Issue 1, 2007, p. 312-329. [Publisher]
  12. Liao S. J., Chwang A. Application of homotopy analysis method in nonlinear oscillations. Journal of Applied Mechanics, Vol. 65, Issue 4, 1998, p. 914-922. [Publisher]
  13. 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. 105-116. [Publisher]
  14. Liao S. J. On the homotopy analysis method for nonlinear problems. Applied Mathematics and Computation, Vol. 147, Issue 2, 2004, p. 499-513. [Publisher]
  15. Park S. H., Kim J. H. Homotopy analysis method for option pricing under stochastic volatility. Applied Mathematics Letters, Vol. 24, Issue 10, 2011, p. 1740-1744. [Publisher]
  16. Liao S. J. The Proposed Homotopy Analysis Technique for the Solution of Nonlinear Problems. Ph.D. Thesis, Shanghai Jiao Tong University, 1992. [Search CrossRef]
  17. Liao S. J. Beyond Perturbation: Introduction to the Homotopy Analysis Method. CRC Press, 2003. [Publisher]
  18. Van Gorder R. A. The variational iteration method is a special case of the homotopy analysis method. Applied Mathematics Letters, Vol. 45, 2015, p. 81-85. [Publisher]
  19. Nayfeh A. H. Perturbation Methods. John Wiley and Sons, 2008. [Search CrossRef]
  20. Marinca V., Herisanu N. Nonlinear Dynamical Systems in Engineering. Springer Berlin Heidelberg, 2012, p. 9-29. [Publisher]