# Iterative Homotopy Harmonic Balance Approach for Determining the Periodic Solution of a Strongly Nonlinear Oscillator

Authors: Huaxiong Chen; Mingkang Ni
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A novel approach about iterative homotopy harmonic balancing is presented to determine the periodic solution for a strongly nonlinear oscillator. This approach does not depend upon the small/large parameter assumption and incorporates the salient features of both methods of the parameter-expansion and the harmonic balance. Importantly, in obtaining the higher-order analytical approximation, all the residual errors are considered in the process of every order approximation to improve the accuracy. With this procedure, the higher-order approximate frequency and corresponding periodic solution can be obtained easily. Comparison of the obtained results with those of the exact solutions shows the high accuracy, simplicity, and efficiency of the approach. The approach can be extended to other nonlinear oscillators in engineering and physics.


---

## Body

## 1. Introduction

Considerable attention has been paid to the study of nonlinear problems not only in all areas of engineering but also in physics and other disciplines, since world is nonlinear essentially and most phenomena are modeled by nonlinear differential equations. However, in general, it is very difficult to get exact solutions for those problems. Thus, to find approximate analytic solutions to these nonlinear problems has been the desire of many researchers for a long time.During the past few decades, there are many analytical methods used to find approximate solutions of various types of nonlinear equations, such as the variational iteration method [1–4], He’s parameter-expansion method [5–7], harmonic balance method [8–12], homotopy perturbation method [13–17], energy balance method [18, 19], and delta-perturbation [20].Among these analytical methods, the iterative homotopy harmonic balancing [21] (denoted as IHHB for abbreviation) is a novel and effective technique for solving some strongly nonlinear oscillators. This approach incorporates the salient features of both methods of the parameter-expansion and the harmonic balance. By constructing a parameter p which is considered as a “small parameter” in the harmonic balance process, the IHHB does not depend upon the large/small parameter assumption and can get high-order analytical approximations easily. Hence, the application of IHHB can be found in various nonlinear problems [21–23].In this paper, we consider the following strongly nonlinear oscillator with nonpolynomial term [6](1)y¨+αy+βy3+γy1/3=0,y0=A,y′0=0,where the over-dot denotes the derivative with respect to the time t and A is the amplitude of the oscillation.Whenα=β=0, γ=1, (1) becomes a truly nonlinear oscillator introduced by Mickens [12] and has been studied by many other investigators [11, 16, 17]. When α=1, β=ϵ, γ=0, (1) reduces to a well-known oscillator namely Duffing equation. By simple analysis [6], we know that (1) has periodic solution when α+βA2+γA-2/3>0.To get higher-order analytical and accurate approximations for the frequency and solution of this oscillator, IHHB solving is set up. More importantly, in obtaining the higher-order analytical approximation, all the residual errors are considered in the process of every order approximation to improve the accuracy. Excellent agreement of the approximate solution and frequency with the exact ones has been demonstrated and discussed. As can be seen, the results presented in this paper reveal that IHHB is effective and convenient for some nonlinear oscillators with nonpolynomial terms.


## 2. Basic Ideas of Iterative Homotopy Harmonic Balancing Approach

Consider the following second-order systems:(2)y¨-fy=0,y0=A,y˙0=0,where the over-dot denotes the derivative with respect to the time t and A is the amplitude of the oscillations. For convenience, we assume f(-y)=-f(y).By introducing a new independent variableτ=ωt, then (2) becomes(3)ω2y′′-fy=0,y0=A,y′0=0,where prime denotes the derivative with respect to τ and ω is an angular frequency to be determined.Let the periodic solutiony(τ) of (3) exist and assume that it can be expressed by such a set of base functions(4)cos⁡2m-1τ∣m=1,2,3,….According to (4), the initial approximate periodic solution satisfying initial conditions should be(5)y0τ=Acos⁡τ,τ=ω0t,where ω0 is an unknown constant to be determined later.Here, we will use the residual to improve the accuracy. Substituting (5) into (3), the initial residual can be obtained:(6)R0τ=ω02y0′′-fy0.IfR0(τ)=0, then y0(τ) happens to be the exact solution. However, in general, such case will not appear for nonlinear problems.According to harmonic balance method, no secular term in (6) requires eliminating contributions proportional to cos⁡τ, through which the unknown constant ω0 can be determined. Then, the zero-order approximation solution is obtained of the form(7)y0τ=Acos⁡τ,τ=ω0t,with initial residual(8)R0τ=ω02y0′′-fy0.In the following, we introduce a bookkeeping parameterp with values in the interval [0,1], denote the y(τ), ω as y(τ,p), ω(p), and then expand y(τ,p) and ω2(p) in a series of p(9)yτ,p=yk-1τ,p+pkykτ+⋯,ω2p=ωk-12p+pkωk+⋯,in which(10)yk-1τ,p=y0τ+py1τ+⋯+pk-1yk-1τ,ωk-12p=ω02+pω1+⋯+pk-1ωk-1.Obviously, whenp=0, y(t,0) is the zeroth-order approximation and when p=1, y(t,1) is the required approximate solution of (3). Generally, (9) provides the higher-order approximation to the exact solution. For example, the first-order analytical approximation turns out to be(11)yτ=y0τ+py1τ,ω2=ω02+pω1.Substituting (11) into (3) and equating the coefficients of the p yield(12)F1τ,ω1,y1τ=ω02y1′′+ω1y0′′-dfdyy0y1.By considering (4) and harmonic balance method, let(13)y1τ=a3,1cos⁡τ-cos⁡3τ.Substituting (13) into (12), we consider the following equation:(14)R1τ=F1τ,ω1,y1τ+R0τ;in this way, the initial residual R0(τ) is introduced into (14) to improve the accuracy.According to harmonic balance method, the right-hand side of (14) should not contain the terms cos⁡τ, cos⁡3τ. Letting their coefficients be zeros yields two linear equations with two unknowns ω1 and a3,1, through which the two unknowns can be solved easily. We now get the first-order approximation(15)y1τ=y0τ+y1τ,ω12=ω02+ω1,τ=ω1t,with the residual(16)R1τ=ω12y1′′-fy1,where y0(τ) and y1(τ) are given by (7) and (13), respectively.With the process going on, the higher-order approximation can be obtained similarly. For example, for thekth-order approximation, we may assume(17)ykτ=∑i=1ka2i+1,kcos⁡τ-cos⁡2i+1τ,k=2,3,….To determine the unknown parametersa2i+1,k(i=2,…,k) and ωk, we substitute (9) into (3) and collect the coefficients of p and finally we get(18)Fkτ,ωk,ykτ=ω02yk′′+∑i=1kωiyk-i′′+1k!dkdpkfyp=0.Substitutingy(k-1)(τ,1) and ω(k-1)2(1) into (3), one yields the following residual:(19)Rk-1τ=ωk-121yk-1′′τ,1-fyk-1τ,1.Let us consider the equation(20)Rkτ=Fkτ,ωk,ykτ+Rk-1τand eliminate the contributions proportional to cos⁡τ,cos⁡3τ,… and cos⁡(2k+1)τ in the right-hand side of (20); we can solve the unknown a3,k,a5,k,…,a2k+1,k and ωk from those equations.Then, thek-order approximations to the frequency and periodic solution can be obtained in the form(21)ykτ=yk-1τ,1+ykτ,ωk2=ωk-121+ωk,where y(k-1)(τ,1) and ω(k-1)2(1) are given by (10).


## 3. Solution Procedure

By introducing a new time scaleτ=ωt, then (1) becomes(22)ω2y′′+αy+βy3+γy1/3=0,y0=A,y′0=0,where prime denotes the derivative with respect τ and ω is an angular frequency to be determined.By following (7), the initial approximation with the initial conditions can be written into(23)y0τ=Acos⁡τ,ω=ω0.Here, it is possible to do the following Fourier series expansion:(24)cos⁡τ1/3=∑n=0∞a2n+1cos⁡2n+1τ,in which(25)a2n+1=4π∫0π/2cos⁡τ1/3cos⁡2n+1τdτ=3Γ7/324/3Γn+5/3Γ2/3-n.Substituting (23) and (24) into (22) yields(26)R0τ=F01cos⁡τ+F03cos⁡3τ+F05cos⁡5τ+⋯,in which(27)F01=αA-ω02A+0.75βA3+1.1596γA1/3,F03=0.25βA3-0.231919γA1/3,F05=0.11596γA1/3.Eliminating contributions proportional tocos⁡τ in (26), that is, solving equation F01=0, gives(28)ω0=α+0.75βA2+1.1596γA-2/3,which is the same as the one by Xu [24]. Till now we get the zeroth-order analytical approximation (23) and the initial residual(29)R0τ=F03cos⁡3τ+F05cos⁡5τ+⋯.To obtain the first-order analytical solution, we substitute (11) and (13) into (12), expand the function (cos⁡τ-cos⁡3τ)(cos⁡τ)-2/3 into a Fourier series(30)cos⁡τ-cos⁡3τcos⁡τ-2/3=∑n=0∞b2n+1cos⁡2n+1τ,in which(31)b2n+1=4π∫0π/2cos⁡τ-cos⁡3τcos⁡τ-2/3cos⁡2n+1τdτ,and finally get(32)F1τ,ω1=F11cos⁡τ+F13cos⁡3τ+F15cos⁡5τ16+⋯,in which(33)F11=0.75a3,1βA2+ω1A-0.695757a3,1γA-2/3,F13=a3,18α+6βA2+9.85656γA-2/3,F15=-0.75a3,1βA2+0.17921a3,1γA-2/3.According to (16), we have(34)R1τ=F1τ,ω1+R0τ.Setting the coefficients ofcos⁡τ, cos⁡3τ to zero in the right-hand side of (34) gives(35)0=0.75a3,1βA2+ω1A-0.695757a3,1γA-2/3,0=0.25βA3-0.231919γA1/3+a3,18α+6βA2+9.85656γA-2/3.Solving (35) yields(36)ω1=-0.0062497112.0006βA8/3-11.1326γ2A2/338.4018αA2/3+28.8013βA8/3+47.3137γ,(37)a3,1=1.11326γA-1.20006βA11/338.4018αA2/3+28.8013βA8/3+47.3137γ.Thus, we obtain the first-order approximation(38)ω1=ω02+ω1,y1t=A+a3,1cos⁡ω1t-a3,1cos⁡3ω1t,with the residual(39)R1τ=0.11596γA1/3-0.75a3,1βA2+0.17921a3,1γA-2/3+⋯,where ω0, ω1, and a3,1 are determined by (28), (36), and (37), respectively.With the procedure going on, similarly, we can get the high-order approximation. For example, the second-order approximation can be obtained as follows:(40)y2t=A+a3,1+a3,2+a5,2cos⁡ω2t-a3,1+a3,2cos⁡3ω2t-a5,2cos⁡5ω2t,(41)ω2=ω02+ω1+ω2,in which(42)a3,2=GL2M,a5,2=HL2M,ω2=FA2/3L2M.The variables L, M, G, H, and F are presented in Appendix.


## 4. Results and Discussion

In order to illustrate and verify the efficiency and correctness of the presented approach for this strongly nonlinear oscillator, we consider some special cases.Case 1.
Ifα=0, β=0, γ=1, (1) reduces to(43)y¨+y1/3=0,y0=A,y′0=0.Then, from (41), the second-order approximation of the frequency can be obtained:(44)ω2=1.1466A-2/3=1.07079A-1/3.The exact frequency [16] is ω=1.070451A-1/3. Therefore, it can be easily proved that the maximal relative error is less than 0.032%.
Hence, from (40), the second-order approximate solution can be expressed as follows:(45)y2t=1.01966Acos⁡ω2t-0.0238744Acos⁡3ω2t+0.00421451Acos⁡5ω2t,which agrees very well with the exact solution [6] as shown in Figure 1.Figure 1
Comparison of the approximate solution with the exact solution;•••• exact solution, — 2nd-order approximate solution.To further illustrate and verify the accuracy of the presented approach in this case, we present the comparison between the approximate and exact frequencies for the second-order approximation by using different methods in Table1. It is clear that, for the second-order approximation, the result obtained in this paper is better than those obtained previously by other authors.Table 1
Comparison of the exact and approximate frequencies (2nd-order) obtained by using different methods.
This paper
HPM [16]
HPM [17]
HBM [11]
HBM [12]
A
1
/
3
ω
2
1.07079
1.06861
1.06991
1.06928
1.06341
% error
0.032
%
0.17
%
0.050
%
0.11
%
0.66
%Case 2.
Ifα=1, β=ϵ, γ=0, (1) reduces to a Duffing equation(46)y¨+y+ϵy3=0,y0=A,y′0=0.
In this case, the first- and second-order approximate periods obtained by (38) and (41) for different values of ϵA2 are shown in Table 2. The relative errors (RE) are defined as RE=|(T-Tex)/Tex|×100. It can be observed that the approximate periods have a good adjustment with the exact ones.Table 2
Comparison of approximate periods with the exact period for Case2.
ϵ
A
2
T
(
1
) (RE)
T
(
2
) (RE)
T
e
x
0.01
6.25976 (0.00000)
6.25976 (0.00000)
6.25976
0.1
6.06066 (0.00000)
6.06066 (0.00000)
6.06066
1
4.76792 (0.00210)
4.76803 (0.00021)
4.76802
10
2.19094 (0.03878)
2.19191 (0.00548)
2.19179
100
0.75814 (2.97178)
0.73633 (0.00951)
0.73626
1000
0.234194 (0.06230)
0.23437 (0.01281)
0.23434To verity results, Figure2 shows the comparisons of the second-order analytical solutions obtained by (40) with the exact ones for A=0.1, 1, 10, and 100 when ϵ=1. It can be seen from this figure that our analytical results are very close to the exact ones for the wide range of initial amplitude in this case.Figure 2
Comparison of the approximate solution with the exact solution;•••• exact solution, — 2nd-order approximate solution.Case 3.
Ifα=1, β=0, γ=1, we can obtain the following nonlinear oscillator:(47)y¨+y+y1/3=0,y0=A,y′0=0.
From Table3 and Figure 3, we can see that the accuracy of the results obtained in this paper is in excellent agreement with exact ones for the wide range of initial amplitude A.Table 3
Comparison of approximate periods with the exact period for Case3.
A
2
T
(
1
) (RE)
T
(
2
) (RE)
T
e
x
0.01
1.24079 (0.10731)
1.23909 (0.02985)
1.23946
0.1
2.49977 (0.08246)
2.49714 (0.02282)
2.49771
1
4.28454 (0.03385)
4.28271 (0.00887)
4.28309
10
5.62191 (0.00481)
5.62157 (0.00125)
5.62164
100
6.12075 (0.00033)
6.12072 (0.00016)
6.12073
1000
6.24708 (0.00016)
6.24707 (0.00000)
6.24707Figure 3
Comparison of the approximate solution with the exact solution;•••• exact solution, — 2nd-order approximate solution.


## 5. Conclusions

An iterative homotopy harmonic balance approach has been presented and applied to deduce the accurate approximations to the angular frequency and periodic solution of a strongly nonlinear oscillator. The high-order analytical approximations of the frequency and solution of this oscillator are obtained. Excellent agreement between approximate results by this approach and the exact ones has been demonstrated and discussed; the discrepancy between the second-order approximate results and exact ones is very very low. And, we can see that the approach considered here is very simple in its principle and has great potential to be applied to other nonlinear oscillators.



---

*Source: 675176-2015-10-28.xml*