- Email: [email protected]
Application of the optimal homotopy asymptotic method for solving a strongly nonlinear oscillatory system
Mathematical and Computer Modelling 58 (2013) 1837–1843
Contents lists available at ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Application of the optimal homotopy asymptotic method for solving a strongly nonlinear oscillatory system A. Golbabai a,∗ , M. Fardi b , K. Sayevand c a
Department of Mathematics, Islamic Azad University, Karaj Branch, Karaj, Iran
b
Department of Mathematics, Islamic Azad University, Boroujen Branch, Boroujen, Iran
c
Department of Mathematics, Faculty of Science, Malayer University, Malayer, Iran
article
abstract
info
Article history: Received 25 April 2011 Received in revised form 23 November 2011 Accepted 5 December 2011 Keywords: Optimal homotopy asymptotic method Homotopy analysis method Homotopy-pade technique
In this paper, the optimal homotopy asymptotic method (OHAM) and the traditional homotopy analysis method (HAM) are used to obtain analytical solution for a strongly nonlinear oscillation. Moreover, the homotopy-pade technique is employed to accelerate the convergence of solution series of traditional HAM. Results show that the second-order approximation by the OHAM is quick convergence and more accurate than the high-order of approximation by the homotopy-pade technique. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Most phenomena in our world are fundamentally nonlinear and formulated by nonlinear equations. Sometimes it is difficult to obtain the accurate solutions of nonlinear equations. There are some analytic approaches for solving nonlinear problems such as the perturbation method [1], the Lyapunov small parameter method [2], the δ -expansion method [3] and the Adomian decomposition method [4] that are well-known. These are widely applied for weakly nonlinear problems. These methods cannot always guarantee the convergence of approximation series. Liao (1992) proposed the homotopy analysis method (HAM) to overcome the restrictions of mentioned traditional techniques [5]. For a given nonlinear equation N (y(t )) = 0,
t ∈ D,
(1)
where N is a nonlinear operator, and y(t ) is an unknown variable, Liao constructed a one-parameter family of equations [5]
(1 − q)L[ϕ(τ ; q) − y0 (t )] + qN [ϕ(τ ; q)] = 0,
t ∈ D, q ∈ [0, 1],
(2)
clearly, we have ϕ(t ; 0) = y0 (t ) and ϕ(t ; 1) = y(t ). Therefore, if the Taylor’s series
ϕ(τ ; q) = y0 (t ) +
∞
yj (t )qj ,
(3)
j =1
converges at q = 1, we have the series solution y(t ) = y0 (t ) +
∞
yj (t ).
j=1
∗
Corresponding author. Tel.: +98 2173225417; fax: +98 2173225417. E-mail addresses: [email protected] (A. Golbabai), [email protected] (M. Fardi), [email protected] (K. Sayevand).
0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.12.027
(4)
1838
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
Liao demonstrated that this early homotopy analysis method cannot always provide us with a convenient way to control the convergence region and rate of approximation series [5–7]. To overcome this restriction, Liao (1997) constructed a two-parameter family of equations [5].
(1 − q)L[ϕ(τ ; q) − y0 (t )] + CqN [ϕ(τ ; q)] = 0,
t ∈ D, q ∈ [0, 1].
(5)
Liao found that the auxiliary parameter C can adjust the convergence region and rate of approximation series [5]. Liao investigated the influence of C on the convergence of solution series by means of C -curves [5], but the C -curves approach cannot give the optimal value of C . Yabushita et al. proposed a method, called ‘‘the optimization method’’ to calculate optimal auxiliary parameter by means of the least square method [8]. Recently, Marinca et al. introduced a new method as the optimal homotopy asymptotic method (OHAM) [9]. They considered such a family of equations
(1 − q)L[ϕ(τ ; q) − y0 (t )] + H (q)N [ϕ(τ ; q)] = 0,
t ∈ D, q ∈ [0, 1],
(6)
where, H (q) is nonzero auxiliary function for q ̸= 0, and the Taylor’s series H (q) = C1 q + C1 q2 + · · · ,
H (0) = 0,
(7)
is converges at q = 1. Obviously, traditional HAM is special case of OHAM when we have H (q) = Cq [10–12]. Like Yabushita et al., Marinca et al. determined the optimal values of the auxiliary parameters by minimizing the square residual error [13]. The OHAM provides us with a convenient way to adjust optimally the convergence of approximation series. In this paper, we use the OHAM and traditional HAM to obtain periodic solutions of the oscillation of a mass attached to a stretched elastic wire and afterwards we discuss about convergence rate of them. 2. Formulation In this section, we consider the oscillation of a mass attached to a stretched elastic wire. The system oscillates between symmetric bounds [−A, A], and its frequency depends on the amplitude A. The oscillation of a mass attached to a stretched elastic wire is as following d2 y dt 2
+y−
λy
1 + y2
= 0,
0 < λ ≤ 1,
y(0) = A,
y′ (0) = 0,
(8)
where y is an unknown real function and λ is known parameter [14]. Now, we write (8) in the following form
(1 + y2 )
d2 y
2y +
dt 2
d2 y
dt 2
+ (1 − λ2 + y2 )y2 = 0,
y(0) = A,
y′ (0) = 0.
(9)
Assume now that the solutions (9) is periodic with the period that T = 2ωπ where ω is frequency of oscillation. Substituting Ω = ω2 , τ = ωt and y(t ) = Y (τ ) in (9) we have
2 d2 Y d Y (1 + Y 2 ) Ω 2 Y 2Y + Ω 2 + (1 − λ2 + Y 2 )Y 3 = 0, dτ dτ
Y (0) = A,
Y ′ (0) = 0.
(10)
Under the initial conditions, Eq. (8) has the exact frequency [15]
ω=
π
1
2
0
Adx
√
A2 (1 − x2 ) − 2λ( 1 + A2 −
√
1 + A 2 x2 )
.
(11)
Remark 1.
ω≈
√
1−λ
for A ≪ 1 and ω ≈ 1 for A ≫ 1.
(12)
3. Application of OHAM To illustrate the basic ideas of the OHAM for differential equations, let Ω0 denote the initial guess of Ω . We choose the auxiliary linear operator as L[ϕ(τ ; q)] = Ω02
d2 ϕ(τ ; q) dτ 2
+ ϕ(τ ; q) .
(13)
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
1839
On the basis (10), we define a nonlinear operator as
N [ϕ(τ ; q), Ω (q)] = (1 + ϕ(τ ; q) ) Ω (q) 2
d2 ϕ(τ ; q) dτ 2
d2 ϕ(τ ; q) ϕ(τ ; q) 2ϕ(τ ; q) + Ω (q) dτ 2
+ (1 − λ2 + ϕ(τ ; q)2 )ϕ(τ ; q)3 .
(14)
By the homotopy technique, we construct a homotopy in a more general form
(1 − q)L[ϕ(τ ; q)] = H (q)N [ϕ(τ ; q), Ω (q)],
q ∈ [0, 1],
(15)
such that
∂ϕ(τ ; q) = 0, ∂τ τ =0
ϕ(0; q) = A
(16)
where ϕ(τ ; q) and Ω (q) are unknown functions and H (q) is an auxiliary function so that H (0) = 0 and H (q) ̸= 0 for q ̸= 0. Obviously, traditional HAM is special case of OHAM when we have H (q) = Cq. For q = 0, we have L[Y0 (τ )] = 0,
Y0 (0) = A,
Y0′ (0) = 0.
(17)
According to (16) we have Y0 (τ ) = A cos(τ ).
(18)
Now we choose the auxiliary function in the form H (q) = pC1 + p2 C2 + · · · ,
(19)
where C1 , C2 , . . . are constants to be determined later. When the parameter q increases from 0 to 1, the solution ϕ(τ ; q) varies from Y0 (τ ) to Y (τ ), so does Ω (q) from Ω0 to Ω . If this continuous variation is smooth enough, the Taylor’s series with respect to q can be constructed for ϕ(τ ; q) and Ω (q).
ϕ(τ ; q) = Y0 (τ ) +
∞
Ym qm ,
Ω (q) = Ω0 +
∞
Ωm qm ,
(20)
m=1
m=1
where Ym (τ ) =
1 ∂ m ϕ(τ ; q) m!
∂ qm
, q=0
Ωm =
1 ∂ m Ω (q) m!
∂ qm
.
(21)
q =0
Assume now that [7,5] N [ϕ(τ ; q), Ω (q)] = N0 (Y0 , Ω0 ) + N1 (Y0 , Ω0 , Y1 , Ω1 )q + · · ·
(22)
where N0 (Y0 , Ω0 ) = N (Y0 , Ω0 ),
(23)
and Ni (Y0 , Ω0 , . . . , Yi , Ωi ) =
∂ i N [ϕ(τ ; q), Ω (q)] |q=0 . i!∂ qi
(24)
In view of homotopy technique, we obtain the following equations. The first and second order equations are given respectively
L[Y1 (τ )] − L[Y0 (τ )] − C1 N0 (Y0 , Ω0 ) = 0, Y1 (0) = 0, Y1′ (0) = 0,
(25)
L[Y2 (τ )] − L[Y1 (τ )] − C1 N1 (Y0 , Ω0 , Y1 , Ω1 ) − C2 N0 (Y0 , Ω0 ) = 0, Y1 (0) = 0, Y1′ (0) = 0.
(26)
and
The general governing equation for Ym (τ ) are given by
m L[Y (τ )] − L[Y (τ )] − C N (Y , Ω , . . . , Y , Ω ) = 0, m m−1 j m−j 0 0 m−j m−j j =1 Y (0) = 0, Y ′ (0) = 0, m ≥ 3. m m
(27)
1840
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
Note that Ym (τ ) and Ωm−1 are all known and we have only Eq. (27) for calculating Ym (τ ). Thus additional algebraic equation is required for determining Ωm−1 . If the term m
Cj Nm−j (Y0 , Ω0 , . . . , Ym−j , Ωm−j ),
(28)
j =1
contains the term cos(τ ) then the solution of linear equation mentioned above involves the so-called secular term τ cos(τ ). Thus coefficient of cos(τ ) must be needed zero. It is to be noted that obviously, free oscillations of a conservation system with odd nonlinearity can be expressed by the base functions
{cos(mωt ) : m = 1, 2, 3, . . .},
(29)
in the form u(τ ) =
∞
ck cos(kτ ),
(30)
k=1
where ck is a coefficient. This provides us with the so-called rule of solution expression. We determine Ωm−1 by the analytic approach mentioned above for m = 1 as
Ω0 = 1 ±
(36 + 30A2 )λ2 . 6 + 5A2
(31)
Since an increase in the amplitude A results in increasing Ω as discussed earlier in formulation section, the negative sign in the equation mentioned above should be selected, i.e.
Ω0 = 1 −
(36 + 30A2 )λ2 . 6 + 5A2
(32)
Also for m = 2, we have
Ω1 =
1 A4 λ2 (−654A2 λ2 + 1100A2 − 636 + 475A4 )
γ
64
C1 ,
where
γ = −150A4 λ2 − 360A2 λ2 − 216λ2 + 450A4 + 125A6 + 540A2 + (36λ2 − 60A2 + 30A2 λ2 − 25A4 − 36) (36 + 30A2 )λ2 . Here, we obtain the constants C1 , C2 , . . . using the least square method. To achieve this, define the vector C⃗M = {C1 , C2 , . . . , CM }.
(33)
In theory at the Mth-order of approximation, we can define the exact square residual error
∆M (C⃗M ) =
M N
D
Yk (τ ),
k=0
M
2 Ωk
dτ .
(34)
k=0
Clearly, the more rapidly ∆M (C⃗M ) decreases to zero, the faster the approximation series converges. So, we can determine the constants Cj as
∂ ∆M = 0, ∂ Cj
j = 1, . . . , M .
(35)
4. Numerical data and results In this section numerical results and discussion about convergence rate of OHAM and traditional HAM is presented. In Tables 1 and 2, we list the required CPU time for calculating the optimal auxiliary parameter of HAM and compared the 10th-order of approximations of Ω with exact solutions (ΩE ). Also, In Tables 3 and 4, we list the required CPU time for calculating the optimal auxiliary parameters of OHAM and compared the different order of approximations of Ω with exact solutions ΩE . As, we see in these tables, it is clear that second-order of approximations by the OHAM is more accurate than 10th-order of approximations by the HAM. In Table 5 we have compared the minimum value of exact square residual error given by 10th-order of approximations (HAM) with the second-order of approximations (OHAM). As we see in these Tables, the minimum values of exact square residual error given by second-order of approximations (OHAM) is less than the minimum values of exact square residual error given by high-order of approximations (HAM).
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
1841
Table 1 The required CPU time for calculating the optimal auxiliary parameter of HAM and comparison of approximations (Ω ) with exact values for (λ, A) = (0.5, 1) and (λ, A) = (0.5, 0.1). M
(λ, A) = (0.5, 1) C
1 2 3 4 5 6 7 8 9 10 11
−2.9174 × 10 −2.2846 × 10−1 −2.1513 × 10−1 −1.3555 × 10−1 −9.3491 × 10−2 −6.9107 × 10−2 −5.3369 × 10−2 −4.2546 × 10−2 −3.4738 × 10−2 −2.8914 × 10−2 −2.4452 × 10−2 −1
(λ, A) = (0.5, 0.1)
(Ω0 +···+ΩM ) ΩE
CPU (s)
C
(Ω0 +···+ΩM ) ΩE
CPU (s)
1.01567 1.01097 1.00623 1.00610 1.00569 1.00577 1.00560 1.00545 1.00532 1.00521 1.00512
2.3 6.3 19.5 33.7 65.1 66.2 143.1 253.9 379.2 640.5 926.6
−27.1593 −15.2916 −9.6916 −6.6545 −4.8383 −3.6709 −2.8782 −2.3160 −1.9033 −1.5915 −1.3503
1.00030 1.00023 1.00019 1.00016 1.00015 1.00013 1.00012 1.00011 1.00011 1.00010 1.00010
5.0 12.5 33.5 55.8 72.7 92.8 126.1 139.3 240.4 330.9 491.4
Table 2 The required CPU time for calculating the optimal auxiliary parameter of HAM and comparison of approximations (Ω ) with exact values for (λ, A) = (0.5, 4) and (λ, A) = (0.5, 0.4). M
(λ, A) = (0.5, 4) C
1 2 3 4 5 6 7 8 9 10 11
−6.6501 × 10 −3.6151 × 10−3 −2.2462 × 10−3 −1.5254 × 10−3 −1.1013 × 10−3 −8.3201 × 10−4 −6.5042 × 10−4 −5.2222 × 10−4 −4.2851 × 10−4 −3.5775 × 10−4 −3.0335 × 10−4 −4
(λ, A) = (0.5, 0.4)
(Ω0 +···+ΩM ) ΩE
CPU (s)
C
(Ω0 +···+ΩM ) ΩE
CPU (s)
1.01742 1.01675 1.01639 1.01618 1.01604 1.01593 1.01586 1.01579 1.01574 1.01571 1.01568
3.3 8.9 17.1 37.4 52.8 93.8 142.4 164.9 233.1 440.8 575.7
−1.7238 −9.9113 × 10−1 −6.3617 × 10−1 −4.4028 × 10−1 −3.2179 × 10−1 −2.4505 × 10−1 −1.9267 × 10−1 −1.5536 × 10−1 −1.2876 × 10−1 −1.0707 × 10−1 −9.0934 × 10−2
1.00422 1.00330 1.00275 1.00241 1.00215 1.00198 1.00185 1.00175 1.00166 1.00159 1.00153
3.7 9.3 12.2 36.2 77.7 90.8 114.9 225.0 271.3 548.9 695.2
Table 3 The required CPU time for calculating the optimal auxiliary parameter of OHAM and comparison the first-order of approximations Ω with exact solutions (ΩE ) for λ = 0.5. A 0.1 0.4 0.7 1 4 7 10
C1
−6.3327 −3.8261 −1.1461 −4.9400 × 10−1 −8.2013 × 10−3 −1.2187 × 10−3 −3.5120 × 10−4
C2
ΩE
−2.6221 −1.8105 −6.5910 × 10−1 −3.4340 × 10−1 −1.5896 × 10−2 −3.6000 × 10−3 −1.3375 × 10−3
5.0186 × 10 5.2726 × 10−1 5.7075 × 10−1 6.1807 × 10−1 8.5143 × 10−1 9.1143 × 10−1 9.3722 × 10−1
−1
(Ω0 +Ω2 ) ΩE
CPU (s)
1.00013 1.00235 1.00700 1.01234 1.01695 1.01133 1.00830
54.7 45.7 58.0 48.4 57.5 52.7 51.2
Table 4 The required CPU time for calculating the optimal auxiliary parameter of OHAM and comparison the second-order of approximations Ω with exact solutions (ΩE ) for λ = 0.5. A
C1
C2
C3
(Ω0 +Ω2 ) ΩE
CPU (s)
0.1 0.4 0.7 1 4 7 10
−6.4477 −3.8734 −6.4290 × 10−1 −1.9400 × 10−1 −5.9210 × 10−4 −2.1762 × 10−4 −1.4577 × 10−4
−1.8062 −1.2243 −6.0130 × 10−1 −4.0060 × 10−1 −2.0713 × 10−2 −5.5754 × 10−3 −4.4003 × 10−4
−2.5651 −1.7923 −3.1970 × 10−1 5.3204 × 10−2 1.6287 × 10−2 2.5710 × 10−2 −9.7053 × 10−4
1.00005 1.00102 1.00477 1.00730 1.01290 1.00734 1.00784
166.7 166.5 167.3 165.5 168.1 169.5 169.7
4.1. Homotopy-pade technique There exist some techniques to improve the convergence rate of a given series by HAM. Among these techniques, the so-called pade technique is widely applied. The so-called homotopy-pade technique was suggested by means of combining
1842
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843 Table 5 Comparison the second-order approximation of OHAM with approximation of homotopy-pade technique for λ = 0.5. Ω[1,1] ΩE Ω[2,2] ΩE Ω[3,3] ΩE Ω[4,4] ΩE Ω[5,5] ΩE Ω[6,6] ΩE Ω[7,7] ΩE ΩOHAM ΩE
A = 0.1
A = 0.4
A=1
A=4
1.00041
1.00593
1.0282
1.01942
1.00013
1.00229
1.01228
1.01794
1.00012
1.00250
1.01327
1.01755
1.00014
1.00265
1.01097
1.01758
1.00013
1.00220
1.01048
1.01750
1.00010
1.00139
1.00948
1.01739
1.00010
1.00122
1.00863
1.01734
1.00005
1.00102
1.00730
1.01290
Table 6 Comparison the minimum value of exact square residual error given by 10th-order of approximations (HAM) with secondorder of approximation (OHAM) for different values of A.
A = 0.1
λ
0.45
0.75
0.95
C1 C2 C3
−75.3524 −24.7635 −37.7306 1.1786 × 10−16 −2.4756 1.4177 × 10−15
−24.5025 −2.4637 −2.1628 2.1694 × 10−15 −5.5764 × 10−1 1.4216 × 10−14
−4.3564 −9.7724 × 10−2 1.2655 × 10−1 6.7768 × 10−15 −9.0732 × 10−2 3.6349 × 10−14
−4.4490 −1.5992 −2.5217 1.3826 × 10−10 −1.6882 × 10−1 1.0286 × 10−9
−1.7079 −2.9823 × 10−1 2.6201 × 10−1 8.6087 × 10−10 −4.0634 × 10−2 1.0601 × 10−8
−5.6324 × 10−1 −1.0364 × 10−1 −3.6065 × 10−2 2.2234 × 10−9 −1.1786 × 10−2 2.2217 × 10−8
−6.7394 × 10−1 −8.1642 × 10−1 −3.5416 × 10−1 7.5505 × 10−8 −8.0451 × 10−2 9.8673 × 10−8
−5.9394 × 10−1 −1.4230 × 10−1 −1.4663 × 10−1 3.1244 × 10−7 −1.6224 × 10−2 1.7771 × 10−6
−1.1378 × 10−1 −2.0008 × 10−1 8.7178 × 10−2 1.0793 × 10−6 −6.8365 × 10−3 4.3473 × 10−6
−2.1622 × 10−1 −5.7148 × 10−1 2.2582 × 10−1 2.1818 × 10−6 −3.7210 × 10−2 1.8221 × 10−6
−1.1647 × 10−1 −2.2758 × 10−1 5.9796 × 10−2 1.4681 × 10−5 −9.5547 × 10−3 3.7126 × 10−5
−1.2588 × 10−1 −4.2276 × 10−2 3.2358 × 10−2 5.6332 × 10−5 −4.6953 × 10−3 1.0451 × 10−4
∆3 (C⃗3 ) C
∆10 (C )
A = 0.4
C1 C2 C3
∆3 (C⃗3 ) C
∆10 (C )
A = 0.7
C1 C2 C3
∆3 (C⃗3 ) C
∆10 (C ) C1 C2 C3
A=1
∆3 (C⃗3 ) C
∆10 (C )
the pade technique with HAM [5]. For a given series M (q) = Ωj qj ,
(36)
j =0
M
the corresponding [k, n] pade approximate is expressed by
WkU,n (q) = (q) = D Wk,n (q) k+n
k
WjU qj
j =0
1+
n
,
(37)
WjD qj
j =1
where WjU
and WjD
are determined by the coefficients Ωj , (j = 0, 1, , n+k). Setting q = 1 provides the [k, n] homotopy-pade
approximation
Ω[k,n]
n+k WkU,n (1) = (1) = Ωj = D = Wk,n (1) k+n j =0
k
WjU
j =0
1+
n j =1
, WjD
(38)
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
1843
which accelerate the convergence rate of solution series of HAM. We have applied the homotopy-pade technique to accelerate the convergence rate of Mth-order approximations of HAM. Remark 2. In Table 6, we compared the approximations of homotopy-pade technique with second-order approximations of OHAM. All of the obtained results show that the second-order approximations of OHAM are more accurate than the approximations of homotopy-pade. 5. Conclusions In this paper, the optimal homotopy asymptotic method (OHAM) and the traditional homotopy analysis method (HAM) were used to obtain analytical solution for a strongly nonlinear oscillation. Moreover, the approximate solution obtained here is valid not only for weakly nonlinear equations, but also for highly nonlinear ones. The convergence and low error is remarkable, which represents that the OHAM has a very high accuracy and it has been observed that OHAM reveals very good agreement with numerical results. According to obtained results, the OHAM does not need to high-order of approximation for solving strongly nonlinear oscillation and conclude that the obtained results of OHAM is more accurate than the approximation of HAM. Acknowledgements The authors are pleased to acknowledge the helpful comments of the referees on this paper. References [1] A.H. Nayfeh, Perturbation Methods, John Wiley & Sons, New York, 2000. [2] A.M. Lyapunov, General Problem on Stability of Motion, Taylor & Francis, London, 1992, (English translation). [3] A.V. Karmishin, A.T. Zhukov, V.G. Kolosov, Methods of dynamics calculation and testing for thin-walled structures, Mashinostroyenie, Moscow, 1990 (in Russian). [4] G. Adomian, Nonlinear stochastic differential equations, J. Math. Anal. Appl. 55 (1976) 441–452. [5] S.J. Liao, Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman and Hall, CRC Press, Boca Raton, 2003. [6] A. Molabahramia, F. Khania, The homotopy analysis method to solve the BurgersHuxley equation, Nonlinear Analysis: RWA 10 (2009) 589–600. [7] M. Turkyilmazoglu, A note on the homotopy analysis method, Appl. Math. Lett. 23 (2010) 1226–1230. [8] K. Yabushita, M. Yamashita, K. Tsuboi, An analytic solution of projectile motion with the quadratic resistance law using the homotopy analysis method, Phys. A: Math. Theor. 40 (2007) 8403–8416. [9] V. Marinca, N. Herisanu, Application of optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer, Int. Commun. Heat Mass Transfer 35 (2008) 710–715. [10] J. Ali, S. Islam, UI. Siraj, Z. Gul, Optimal homotopy asymptotic method for multipoint boundary value problems, Comput. Math. Appl. (2010) 2000–2006. [11] V. Marinca, N. Herisanu, An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate, Appl. Math. Lett. (2009) 245–251. [12] M. Idrees, S. Islam, H. Sirajul, I. Sirajul, Application of optimal homotopy asymptotic method to squeezing flow, Comput. Math. Appl. 59 (11) (2010) 3858–3866. [13] V. Marinca, N. Herisanu, et al., An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate, Appl. Math. Lett. 22 (2009) 245–251. [14] A. Beelendez, A. Hernandez, T. Belendez, M.L. Alvarez, S. Gallego, M. Ortuno, C. Neipp, Application of the harmonic balance method to a nonlinear oscillator typified by a mass attached to stretched wire, J. Sound Vib. 302 (2007) 1018–1029. [15] W.P. Sun, B.S. Wu, C.W. Lim, Approximate analytical solution for oscillation of a mass attached to stretched elastic wire, J. Sound Vib. 300 (2007) 1042–1047.
Contents lists available at ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Application of the optimal homotopy asymptotic method for solving a strongly nonlinear oscillatory system A. Golbabai a,∗ , M. Fardi b , K. Sayevand c a
Department of Mathematics, Islamic Azad University, Karaj Branch, Karaj, Iran
b
Department of Mathematics, Islamic Azad University, Boroujen Branch, Boroujen, Iran
c
Department of Mathematics, Faculty of Science, Malayer University, Malayer, Iran
article
abstract
info
Article history: Received 25 April 2011 Received in revised form 23 November 2011 Accepted 5 December 2011 Keywords: Optimal homotopy asymptotic method Homotopy analysis method Homotopy-pade technique
In this paper, the optimal homotopy asymptotic method (OHAM) and the traditional homotopy analysis method (HAM) are used to obtain analytical solution for a strongly nonlinear oscillation. Moreover, the homotopy-pade technique is employed to accelerate the convergence of solution series of traditional HAM. Results show that the second-order approximation by the OHAM is quick convergence and more accurate than the high-order of approximation by the homotopy-pade technique. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Most phenomena in our world are fundamentally nonlinear and formulated by nonlinear equations. Sometimes it is difficult to obtain the accurate solutions of nonlinear equations. There are some analytic approaches for solving nonlinear problems such as the perturbation method [1], the Lyapunov small parameter method [2], the δ -expansion method [3] and the Adomian decomposition method [4] that are well-known. These are widely applied for weakly nonlinear problems. These methods cannot always guarantee the convergence of approximation series. Liao (1992) proposed the homotopy analysis method (HAM) to overcome the restrictions of mentioned traditional techniques [5]. For a given nonlinear equation N (y(t )) = 0,
t ∈ D,
(1)
where N is a nonlinear operator, and y(t ) is an unknown variable, Liao constructed a one-parameter family of equations [5]
(1 − q)L[ϕ(τ ; q) − y0 (t )] + qN [ϕ(τ ; q)] = 0,
t ∈ D, q ∈ [0, 1],
(2)
clearly, we have ϕ(t ; 0) = y0 (t ) and ϕ(t ; 1) = y(t ). Therefore, if the Taylor’s series
ϕ(τ ; q) = y0 (t ) +
∞
yj (t )qj ,
(3)
j =1
converges at q = 1, we have the series solution y(t ) = y0 (t ) +
∞
yj (t ).
j=1
∗
Corresponding author. Tel.: +98 2173225417; fax: +98 2173225417. E-mail addresses: [email protected] (A. Golbabai), [email protected] (M. Fardi), [email protected] (K. Sayevand).
0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.12.027
(4)
1838
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
Liao demonstrated that this early homotopy analysis method cannot always provide us with a convenient way to control the convergence region and rate of approximation series [5–7]. To overcome this restriction, Liao (1997) constructed a two-parameter family of equations [5].
(1 − q)L[ϕ(τ ; q) − y0 (t )] + CqN [ϕ(τ ; q)] = 0,
t ∈ D, q ∈ [0, 1].
(5)
Liao found that the auxiliary parameter C can adjust the convergence region and rate of approximation series [5]. Liao investigated the influence of C on the convergence of solution series by means of C -curves [5], but the C -curves approach cannot give the optimal value of C . Yabushita et al. proposed a method, called ‘‘the optimization method’’ to calculate optimal auxiliary parameter by means of the least square method [8]. Recently, Marinca et al. introduced a new method as the optimal homotopy asymptotic method (OHAM) [9]. They considered such a family of equations
(1 − q)L[ϕ(τ ; q) − y0 (t )] + H (q)N [ϕ(τ ; q)] = 0,
t ∈ D, q ∈ [0, 1],
(6)
where, H (q) is nonzero auxiliary function for q ̸= 0, and the Taylor’s series H (q) = C1 q + C1 q2 + · · · ,
H (0) = 0,
(7)
is converges at q = 1. Obviously, traditional HAM is special case of OHAM when we have H (q) = Cq [10–12]. Like Yabushita et al., Marinca et al. determined the optimal values of the auxiliary parameters by minimizing the square residual error [13]. The OHAM provides us with a convenient way to adjust optimally the convergence of approximation series. In this paper, we use the OHAM and traditional HAM to obtain periodic solutions of the oscillation of a mass attached to a stretched elastic wire and afterwards we discuss about convergence rate of them. 2. Formulation In this section, we consider the oscillation of a mass attached to a stretched elastic wire. The system oscillates between symmetric bounds [−A, A], and its frequency depends on the amplitude A. The oscillation of a mass attached to a stretched elastic wire is as following d2 y dt 2
+y−
λy
1 + y2
= 0,
0 < λ ≤ 1,
y(0) = A,
y′ (0) = 0,
(8)
where y is an unknown real function and λ is known parameter [14]. Now, we write (8) in the following form
(1 + y2 )
d2 y
2y +
dt 2
d2 y
dt 2
+ (1 − λ2 + y2 )y2 = 0,
y(0) = A,
y′ (0) = 0.
(9)
Assume now that the solutions (9) is periodic with the period that T = 2ωπ where ω is frequency of oscillation. Substituting Ω = ω2 , τ = ωt and y(t ) = Y (τ ) in (9) we have
2 d2 Y d Y (1 + Y 2 ) Ω 2 Y 2Y + Ω 2 + (1 − λ2 + Y 2 )Y 3 = 0, dτ dτ
Y (0) = A,
Y ′ (0) = 0.
(10)
Under the initial conditions, Eq. (8) has the exact frequency [15]
ω=
π
1
2
0
Adx
√
A2 (1 − x2 ) − 2λ( 1 + A2 −
√
1 + A 2 x2 )
.
(11)
Remark 1.
ω≈
√
1−λ
for A ≪ 1 and ω ≈ 1 for A ≫ 1.
(12)
3. Application of OHAM To illustrate the basic ideas of the OHAM for differential equations, let Ω0 denote the initial guess of Ω . We choose the auxiliary linear operator as L[ϕ(τ ; q)] = Ω02
d2 ϕ(τ ; q) dτ 2
+ ϕ(τ ; q) .
(13)
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
1839
On the basis (10), we define a nonlinear operator as
N [ϕ(τ ; q), Ω (q)] = (1 + ϕ(τ ; q) ) Ω (q) 2
d2 ϕ(τ ; q) dτ 2
d2 ϕ(τ ; q) ϕ(τ ; q) 2ϕ(τ ; q) + Ω (q) dτ 2
+ (1 − λ2 + ϕ(τ ; q)2 )ϕ(τ ; q)3 .
(14)
By the homotopy technique, we construct a homotopy in a more general form
(1 − q)L[ϕ(τ ; q)] = H (q)N [ϕ(τ ; q), Ω (q)],
q ∈ [0, 1],
(15)
such that
∂ϕ(τ ; q) = 0, ∂τ τ =0
ϕ(0; q) = A
(16)
where ϕ(τ ; q) and Ω (q) are unknown functions and H (q) is an auxiliary function so that H (0) = 0 and H (q) ̸= 0 for q ̸= 0. Obviously, traditional HAM is special case of OHAM when we have H (q) = Cq. For q = 0, we have L[Y0 (τ )] = 0,
Y0 (0) = A,
Y0′ (0) = 0.
(17)
According to (16) we have Y0 (τ ) = A cos(τ ).
(18)
Now we choose the auxiliary function in the form H (q) = pC1 + p2 C2 + · · · ,
(19)
where C1 , C2 , . . . are constants to be determined later. When the parameter q increases from 0 to 1, the solution ϕ(τ ; q) varies from Y0 (τ ) to Y (τ ), so does Ω (q) from Ω0 to Ω . If this continuous variation is smooth enough, the Taylor’s series with respect to q can be constructed for ϕ(τ ; q) and Ω (q).
ϕ(τ ; q) = Y0 (τ ) +
∞
Ym qm ,
Ω (q) = Ω0 +
∞
Ωm qm ,
(20)
m=1
m=1
where Ym (τ ) =
1 ∂ m ϕ(τ ; q) m!
∂ qm
, q=0
Ωm =
1 ∂ m Ω (q) m!
∂ qm
.
(21)
q =0
Assume now that [7,5] N [ϕ(τ ; q), Ω (q)] = N0 (Y0 , Ω0 ) + N1 (Y0 , Ω0 , Y1 , Ω1 )q + · · ·
(22)
where N0 (Y0 , Ω0 ) = N (Y0 , Ω0 ),
(23)
and Ni (Y0 , Ω0 , . . . , Yi , Ωi ) =
∂ i N [ϕ(τ ; q), Ω (q)] |q=0 . i!∂ qi
(24)
In view of homotopy technique, we obtain the following equations. The first and second order equations are given respectively
L[Y1 (τ )] − L[Y0 (τ )] − C1 N0 (Y0 , Ω0 ) = 0, Y1 (0) = 0, Y1′ (0) = 0,
(25)
L[Y2 (τ )] − L[Y1 (τ )] − C1 N1 (Y0 , Ω0 , Y1 , Ω1 ) − C2 N0 (Y0 , Ω0 ) = 0, Y1 (0) = 0, Y1′ (0) = 0.
(26)
and
The general governing equation for Ym (τ ) are given by
m L[Y (τ )] − L[Y (τ )] − C N (Y , Ω , . . . , Y , Ω ) = 0, m m−1 j m−j 0 0 m−j m−j j =1 Y (0) = 0, Y ′ (0) = 0, m ≥ 3. m m
(27)
1840
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
Note that Ym (τ ) and Ωm−1 are all known and we have only Eq. (27) for calculating Ym (τ ). Thus additional algebraic equation is required for determining Ωm−1 . If the term m
Cj Nm−j (Y0 , Ω0 , . . . , Ym−j , Ωm−j ),
(28)
j =1
contains the term cos(τ ) then the solution of linear equation mentioned above involves the so-called secular term τ cos(τ ). Thus coefficient of cos(τ ) must be needed zero. It is to be noted that obviously, free oscillations of a conservation system with odd nonlinearity can be expressed by the base functions
{cos(mωt ) : m = 1, 2, 3, . . .},
(29)
in the form u(τ ) =
∞
ck cos(kτ ),
(30)
k=1
where ck is a coefficient. This provides us with the so-called rule of solution expression. We determine Ωm−1 by the analytic approach mentioned above for m = 1 as
Ω0 = 1 ±
(36 + 30A2 )λ2 . 6 + 5A2
(31)
Since an increase in the amplitude A results in increasing Ω as discussed earlier in formulation section, the negative sign in the equation mentioned above should be selected, i.e.
Ω0 = 1 −
(36 + 30A2 )λ2 . 6 + 5A2
(32)
Also for m = 2, we have
Ω1 =
1 A4 λ2 (−654A2 λ2 + 1100A2 − 636 + 475A4 )
γ
64
C1 ,
where
γ = −150A4 λ2 − 360A2 λ2 − 216λ2 + 450A4 + 125A6 + 540A2 + (36λ2 − 60A2 + 30A2 λ2 − 25A4 − 36) (36 + 30A2 )λ2 . Here, we obtain the constants C1 , C2 , . . . using the least square method. To achieve this, define the vector C⃗M = {C1 , C2 , . . . , CM }.
(33)
In theory at the Mth-order of approximation, we can define the exact square residual error
∆M (C⃗M ) =
M N
D
Yk (τ ),
k=0
M
2 Ωk
dτ .
(34)
k=0
Clearly, the more rapidly ∆M (C⃗M ) decreases to zero, the faster the approximation series converges. So, we can determine the constants Cj as
∂ ∆M = 0, ∂ Cj
j = 1, . . . , M .
(35)
4. Numerical data and results In this section numerical results and discussion about convergence rate of OHAM and traditional HAM is presented. In Tables 1 and 2, we list the required CPU time for calculating the optimal auxiliary parameter of HAM and compared the 10th-order of approximations of Ω with exact solutions (ΩE ). Also, In Tables 3 and 4, we list the required CPU time for calculating the optimal auxiliary parameters of OHAM and compared the different order of approximations of Ω with exact solutions ΩE . As, we see in these tables, it is clear that second-order of approximations by the OHAM is more accurate than 10th-order of approximations by the HAM. In Table 5 we have compared the minimum value of exact square residual error given by 10th-order of approximations (HAM) with the second-order of approximations (OHAM). As we see in these Tables, the minimum values of exact square residual error given by second-order of approximations (OHAM) is less than the minimum values of exact square residual error given by high-order of approximations (HAM).
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
1841
Table 1 The required CPU time for calculating the optimal auxiliary parameter of HAM and comparison of approximations (Ω ) with exact values for (λ, A) = (0.5, 1) and (λ, A) = (0.5, 0.1). M
(λ, A) = (0.5, 1) C
1 2 3 4 5 6 7 8 9 10 11
−2.9174 × 10 −2.2846 × 10−1 −2.1513 × 10−1 −1.3555 × 10−1 −9.3491 × 10−2 −6.9107 × 10−2 −5.3369 × 10−2 −4.2546 × 10−2 −3.4738 × 10−2 −2.8914 × 10−2 −2.4452 × 10−2 −1
(λ, A) = (0.5, 0.1)
(Ω0 +···+ΩM ) ΩE
CPU (s)
C
(Ω0 +···+ΩM ) ΩE
CPU (s)
1.01567 1.01097 1.00623 1.00610 1.00569 1.00577 1.00560 1.00545 1.00532 1.00521 1.00512
2.3 6.3 19.5 33.7 65.1 66.2 143.1 253.9 379.2 640.5 926.6
−27.1593 −15.2916 −9.6916 −6.6545 −4.8383 −3.6709 −2.8782 −2.3160 −1.9033 −1.5915 −1.3503
1.00030 1.00023 1.00019 1.00016 1.00015 1.00013 1.00012 1.00011 1.00011 1.00010 1.00010
5.0 12.5 33.5 55.8 72.7 92.8 126.1 139.3 240.4 330.9 491.4
Table 2 The required CPU time for calculating the optimal auxiliary parameter of HAM and comparison of approximations (Ω ) with exact values for (λ, A) = (0.5, 4) and (λ, A) = (0.5, 0.4). M
(λ, A) = (0.5, 4) C
1 2 3 4 5 6 7 8 9 10 11
−6.6501 × 10 −3.6151 × 10−3 −2.2462 × 10−3 −1.5254 × 10−3 −1.1013 × 10−3 −8.3201 × 10−4 −6.5042 × 10−4 −5.2222 × 10−4 −4.2851 × 10−4 −3.5775 × 10−4 −3.0335 × 10−4 −4
(λ, A) = (0.5, 0.4)
(Ω0 +···+ΩM ) ΩE
CPU (s)
C
(Ω0 +···+ΩM ) ΩE
CPU (s)
1.01742 1.01675 1.01639 1.01618 1.01604 1.01593 1.01586 1.01579 1.01574 1.01571 1.01568
3.3 8.9 17.1 37.4 52.8 93.8 142.4 164.9 233.1 440.8 575.7
−1.7238 −9.9113 × 10−1 −6.3617 × 10−1 −4.4028 × 10−1 −3.2179 × 10−1 −2.4505 × 10−1 −1.9267 × 10−1 −1.5536 × 10−1 −1.2876 × 10−1 −1.0707 × 10−1 −9.0934 × 10−2
1.00422 1.00330 1.00275 1.00241 1.00215 1.00198 1.00185 1.00175 1.00166 1.00159 1.00153
3.7 9.3 12.2 36.2 77.7 90.8 114.9 225.0 271.3 548.9 695.2
Table 3 The required CPU time for calculating the optimal auxiliary parameter of OHAM and comparison the first-order of approximations Ω with exact solutions (ΩE ) for λ = 0.5. A 0.1 0.4 0.7 1 4 7 10
C1
−6.3327 −3.8261 −1.1461 −4.9400 × 10−1 −8.2013 × 10−3 −1.2187 × 10−3 −3.5120 × 10−4
C2
ΩE
−2.6221 −1.8105 −6.5910 × 10−1 −3.4340 × 10−1 −1.5896 × 10−2 −3.6000 × 10−3 −1.3375 × 10−3
5.0186 × 10 5.2726 × 10−1 5.7075 × 10−1 6.1807 × 10−1 8.5143 × 10−1 9.1143 × 10−1 9.3722 × 10−1
−1
(Ω0 +Ω2 ) ΩE
CPU (s)
1.00013 1.00235 1.00700 1.01234 1.01695 1.01133 1.00830
54.7 45.7 58.0 48.4 57.5 52.7 51.2
Table 4 The required CPU time for calculating the optimal auxiliary parameter of OHAM and comparison the second-order of approximations Ω with exact solutions (ΩE ) for λ = 0.5. A
C1
C2
C3
(Ω0 +Ω2 ) ΩE
CPU (s)
0.1 0.4 0.7 1 4 7 10
−6.4477 −3.8734 −6.4290 × 10−1 −1.9400 × 10−1 −5.9210 × 10−4 −2.1762 × 10−4 −1.4577 × 10−4
−1.8062 −1.2243 −6.0130 × 10−1 −4.0060 × 10−1 −2.0713 × 10−2 −5.5754 × 10−3 −4.4003 × 10−4
−2.5651 −1.7923 −3.1970 × 10−1 5.3204 × 10−2 1.6287 × 10−2 2.5710 × 10−2 −9.7053 × 10−4
1.00005 1.00102 1.00477 1.00730 1.01290 1.00734 1.00784
166.7 166.5 167.3 165.5 168.1 169.5 169.7
4.1. Homotopy-pade technique There exist some techniques to improve the convergence rate of a given series by HAM. Among these techniques, the so-called pade technique is widely applied. The so-called homotopy-pade technique was suggested by means of combining
1842
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843 Table 5 Comparison the second-order approximation of OHAM with approximation of homotopy-pade technique for λ = 0.5. Ω[1,1] ΩE Ω[2,2] ΩE Ω[3,3] ΩE Ω[4,4] ΩE Ω[5,5] ΩE Ω[6,6] ΩE Ω[7,7] ΩE ΩOHAM ΩE
A = 0.1
A = 0.4
A=1
A=4
1.00041
1.00593
1.0282
1.01942
1.00013
1.00229
1.01228
1.01794
1.00012
1.00250
1.01327
1.01755
1.00014
1.00265
1.01097
1.01758
1.00013
1.00220
1.01048
1.01750
1.00010
1.00139
1.00948
1.01739
1.00010
1.00122
1.00863
1.01734
1.00005
1.00102
1.00730
1.01290
Table 6 Comparison the minimum value of exact square residual error given by 10th-order of approximations (HAM) with secondorder of approximation (OHAM) for different values of A.
A = 0.1
λ
0.45
0.75
0.95
C1 C2 C3
−75.3524 −24.7635 −37.7306 1.1786 × 10−16 −2.4756 1.4177 × 10−15
−24.5025 −2.4637 −2.1628 2.1694 × 10−15 −5.5764 × 10−1 1.4216 × 10−14
−4.3564 −9.7724 × 10−2 1.2655 × 10−1 6.7768 × 10−15 −9.0732 × 10−2 3.6349 × 10−14
−4.4490 −1.5992 −2.5217 1.3826 × 10−10 −1.6882 × 10−1 1.0286 × 10−9
−1.7079 −2.9823 × 10−1 2.6201 × 10−1 8.6087 × 10−10 −4.0634 × 10−2 1.0601 × 10−8
−5.6324 × 10−1 −1.0364 × 10−1 −3.6065 × 10−2 2.2234 × 10−9 −1.1786 × 10−2 2.2217 × 10−8
−6.7394 × 10−1 −8.1642 × 10−1 −3.5416 × 10−1 7.5505 × 10−8 −8.0451 × 10−2 9.8673 × 10−8
−5.9394 × 10−1 −1.4230 × 10−1 −1.4663 × 10−1 3.1244 × 10−7 −1.6224 × 10−2 1.7771 × 10−6
−1.1378 × 10−1 −2.0008 × 10−1 8.7178 × 10−2 1.0793 × 10−6 −6.8365 × 10−3 4.3473 × 10−6
−2.1622 × 10−1 −5.7148 × 10−1 2.2582 × 10−1 2.1818 × 10−6 −3.7210 × 10−2 1.8221 × 10−6
−1.1647 × 10−1 −2.2758 × 10−1 5.9796 × 10−2 1.4681 × 10−5 −9.5547 × 10−3 3.7126 × 10−5
−1.2588 × 10−1 −4.2276 × 10−2 3.2358 × 10−2 5.6332 × 10−5 −4.6953 × 10−3 1.0451 × 10−4
∆3 (C⃗3 ) C
∆10 (C )
A = 0.4
C1 C2 C3
∆3 (C⃗3 ) C
∆10 (C )
A = 0.7
C1 C2 C3
∆3 (C⃗3 ) C
∆10 (C ) C1 C2 C3
A=1
∆3 (C⃗3 ) C
∆10 (C )
the pade technique with HAM [5]. For a given series M (q) = Ωj qj ,
(36)
j =0
M
the corresponding [k, n] pade approximate is expressed by
WkU,n (q) = (q) = D Wk,n (q) k+n
k
WjU qj
j =0
1+
n
,
(37)
WjD qj
j =1
where WjU
and WjD
are determined by the coefficients Ωj , (j = 0, 1, , n+k). Setting q = 1 provides the [k, n] homotopy-pade
approximation
Ω[k,n]
n+k WkU,n (1) = (1) = Ωj = D = Wk,n (1) k+n j =0
k
WjU
j =0
1+
n j =1
, WjD
(38)
A. Golbabai et al. / Mathematical and Computer Modelling 58 (2013) 1837–1843
1843
which accelerate the convergence rate of solution series of HAM. We have applied the homotopy-pade technique to accelerate the convergence rate of Mth-order approximations of HAM. Remark 2. In Table 6, we compared the approximations of homotopy-pade technique with second-order approximations of OHAM. All of the obtained results show that the second-order approximations of OHAM are more accurate than the approximations of homotopy-pade. 5. Conclusions In this paper, the optimal homotopy asymptotic method (OHAM) and the traditional homotopy analysis method (HAM) were used to obtain analytical solution for a strongly nonlinear oscillation. Moreover, the approximate solution obtained here is valid not only for weakly nonlinear equations, but also for highly nonlinear ones. The convergence and low error is remarkable, which represents that the OHAM has a very high accuracy and it has been observed that OHAM reveals very good agreement with numerical results. According to obtained results, the OHAM does not need to high-order of approximation for solving strongly nonlinear oscillation and conclude that the obtained results of OHAM is more accurate than the approximation of HAM. Acknowledgements The authors are pleased to acknowledge the helpful comments of the referees on this paper. References [1] A.H. Nayfeh, Perturbation Methods, John Wiley & Sons, New York, 2000. [2] A.M. Lyapunov, General Problem on Stability of Motion, Taylor & Francis, London, 1992, (English translation). [3] A.V. Karmishin, A.T. Zhukov, V.G. Kolosov, Methods of dynamics calculation and testing for thin-walled structures, Mashinostroyenie, Moscow, 1990 (in Russian). [4] G. Adomian, Nonlinear stochastic differential equations, J. Math. Anal. Appl. 55 (1976) 441–452. [5] S.J. Liao, Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman and Hall, CRC Press, Boca Raton, 2003. [6] A. Molabahramia, F. Khania, The homotopy analysis method to solve the BurgersHuxley equation, Nonlinear Analysis: RWA 10 (2009) 589–600. [7] M. Turkyilmazoglu, A note on the homotopy analysis method, Appl. Math. Lett. 23 (2010) 1226–1230. [8] K. Yabushita, M. Yamashita, K. Tsuboi, An analytic solution of projectile motion with the quadratic resistance law using the homotopy analysis method, Phys. A: Math. Theor. 40 (2007) 8403–8416. [9] V. Marinca, N. Herisanu, Application of optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer, Int. Commun. Heat Mass Transfer 35 (2008) 710–715. [10] J. Ali, S. Islam, UI. Siraj, Z. Gul, Optimal homotopy asymptotic method for multipoint boundary value problems, Comput. Math. Appl. (2010) 2000–2006. [11] V. Marinca, N. Herisanu, An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate, Appl. Math. Lett. (2009) 245–251. [12] M. Idrees, S. Islam, H. Sirajul, I. Sirajul, Application of optimal homotopy asymptotic method to squeezing flow, Comput. Math. Appl. 59 (11) (2010) 3858–3866. [13] V. Marinca, N. Herisanu, et al., An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate, Appl. Math. Lett. 22 (2009) 245–251. [14] A. Beelendez, A. Hernandez, T. Belendez, M.L. Alvarez, S. Gallego, M. Ortuno, C. Neipp, Application of the harmonic balance method to a nonlinear oscillator typified by a mass attached to stretched wire, J. Sound Vib. 302 (2007) 1018–1029. [15] W.P. Sun, B.S. Wu, C.W. Lim, Approximate analytical solution for oscillation of a mass attached to stretched elastic wire, J. Sound Vib. 300 (2007) 1042–1047.