# Comparative Vibration Analysis of  a Parametrically Nonlinear Excited Oscillator  Using HPM and Numerical Method

Authors: I. Khatami; M. H. Pashai; N. Tolou
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

The objective of this paper is to present an analytical investigation to analyze the vibration of parametrically excited oscillator with strong cubic negative nonlinearity based on Mathieu-Duffing equation. The analytic investigation was conducted by using He's homotopy-perturbation method (HPM). In order to obtain the analytical solution of Mathieu-Duffing equation, homotopy-perturbation method has been utilized. The Runge-Kutta's (RK) algorithm was used to solve the governing equation via numerical solution. Finally, to demonstrate the validity of the proposed method, the response of the oscillator, which was obtained from approximate solution, has been shown graphically and compared with that of numerical solution. Afterward, the effects of variation of the parameters on the accuracy of the homotopy- perturbation method were studied.


---

## Body

## 1. Introduction

Parametrically excited  systems are widely
spread in many branches of physics and engineering. It is important to
investigate the dynamic behavior of these systems.Vibrations of such
systems occur in wide range of mechanics, due to time-varying loads, especially periodic ones.
These vibrations appear in columns made of nonlinear elastic material [1],
beams with a harmonically variable length [2], beams with harmonic motion of
their support [3], floating offshore structures [4], parametrically excited
pendulums [5], cables being towed by a submarine [6, 7], and so forth.
Parametric excitations occur in electrostatically driven
microelectromechanical oscillators [8], which are produced by fluctuating
voltages applied across comb drives. In practical engineering situations, the
properties of parametric oscillations are widely used, for example, in the
radio, the computer, and laser engineering, in vibrant machines with special
design [9], Paul trap mass spectrometers [10], as well as in a simulator for proving the
equivalence of inertia and passive gravitational mass [11]. Parametric
resonance has been well established in many areas of science, including the
stability of ships, the forced motion of a swing, and Faraday surface wave
patterns  on water. The highly sensitive mass sensor is studied as an in-plane
parametrically resonant oscillator [12]. The simplest mathematical model of the
system with a parametric periodic load is usually a linear Mathieu differential
equation. Due to the nonlinear properties of a real system, nonlinear terms are
added to the equation by Younesian et al. [13]. They have added a cubic
nonlinearity term (x3) to Mathieu-Duffing equation
using averaging method [13, 14]. It is obvious that for such nonlinear equations, there
is no precise analytical solution; consequently these nonlinear equations need
to be solved using other methods. In recent decades, numerical methods have
been well established for analyzing the nonlinear equations such as Mathieu-differential
equation. Most scientists believe that the combination of numerical and
semianalytical methods can give rise to useful results. The two common methods employed to solve the Mathieu-differential
nonlinear equation are based on the numerical integration by Jazar [15] and
perturbation by Cveticanin and Kovacic [16]. Perturbation method, which is much more developed,
is one of the well-known methods that were studied by some researchers [17–19] to solve the
linear and nonlinear equations. Actually, these scientists had paid more
attention to the mathematical aspects. Since there are some limitations in
using the common perturbation method together with the fact that this method is
based on the existence of small parameter, developing the method for different
applications is very difficult. Therefore, a number of semianalytical methods,
including the variational iteration method [20, 21] and the homotopy-perturbation method (HPM),
[22–24] have recently been developed by using new techniques to eliminate the
small parameter. Other authors have made much contribution to the applications
of these methods [25–37].In this paper, we extend the work of
Cveticanin and Kovacic [16] by applying HPM to analyze the vibrations of parametrically
excited oscillator with strong cubic negative nonlinearity. At the first, a brief
review of the method will be given and then the HPM will be applied to the Mathieu-Duffing equation
(1.1). Next,
Runge-Kutta’s (RK) algorithm will be introduced to solve the governing equation.
Finally, a numerical example is given to demonstrate the validity of the
proposed method (HPM) and the effect of the variation parameters on the accuracy
of the HPM is investigated.The governing equation of Mathieu-Duffing system which
is considered in this study is described by the following high-order nonlinear
differential equation:(1.1)x¨+[δ+2εcos⁡(2t)]x−ϕx3=0,where dots indicate differentiation with respect to the
time (t), ε≪1 is a small parameter, φ is the parameter of nonlinearity, and δ is the transient curve and can
be defined as [16](1.2)δ=ϕx02(1−2ε2+ϕx02).The initial condition considered in this study is defined by x = 0.1, x˙ = 0 [16].This
mathematical model corresponds to the parametrically excited oscillator with a
softening spring.


## 2. Analytical Solution

The
analytical solution of such high-order nonlinear differential equation as (1.1)
is usually difficult to find. For these equations, much attention was paid to
develop various simple numerical and analytical methods such as those based on
the numerical integration [15] and perturbation [16]. In this paper, an
approximate analytical method is developed to solve the problem. However, as mentioned
before, for this case, some analytical solutions are well established.
Application of the proposed semianalytical method for achieving the
nonlinearity of parametrically excited vibrations of an oscillator is also discussed
in this paper.


### 2.1. Homotopy-Perturbation Method (HPM)

#### 2.1.1. Background

In this paper, we apply HPM [22–37] for the
solution of the Mathieu-Duffing
equation (1.1). In order to demonstrate how this method works, let us consider
the following nonlinear differential equation:(2.1)A(u)−f(r)=0,r∈Ω,subject to the
boundary  conditions  of(2.2)B(u,∂u∂r)=0,r∈Γ,where
A is a general differential operator, B is a boundary operator, f(r)
is a known analytic functional, and Γ is the boundary of the domain Ω.The
operatorA can generally be divided into two parts L and
N, where L is linear, whereas N is nonlinear. Therefore, (2.1)
can be rewritten as(2.3)L(u)+N(u)−f(r)=0.Homotopy-perturbation
structure is shown to be as the following equation [22–37]:(2.4)H(v,p)=(1−p)[L(v)−L(u0)]+p[L(v)+N(v)−f(r)]=0,where(2.5)v(r,p):Ω×[0,1]→R.For
p = 0 and p = 1, (2.4) reduces to the following equations,
respectively:(2.6)H(v,0)=L(v)−L(u0)=0,H(v,1)=A(v)−f(r)=0,where p∈[0,1] is an embedding parameter and u0 is the first approximation that satisfied the
boundary condition.The
process of changes inp from zero to unity is that of v(r,p) changing form u0 to u(r).
We consider v, as the following [22–37]:(2.7)v=v0+pv1+p2v2+⋯,and
the best approximation for solution is [22–37](2.8)u=lim⁡p→1v=v0+v1+v2+⋯.


#### 2.1.2. Implementation of HPM

Now we apply homotopy-perturbation (2.4)
to (1.1). We construct a homotopy in the following form:(2.9)(1−p)[d2x(t)dt2]+p{−[δ+2εcos⁡(2t)]x(t)+ϕx(t)3}=0.Obviously,
(2.9) becomes linear if p = 0, when p = 1, it becomes the original nonlinear
one. So the variation of p from zero to unity makes the equation to change to nonlinear
from linear one.According
to HPM, we assume that the solution of (2.9) can be expressed in a series of p(2.10)v(r)=v0(r)+pv1(r)+p2v2(r)+⋯.Substituting v(r) from (2.10)
into (2.9), after some simplification and substitution and rearranging based on
powers of p-terms, we have(2.11)p0:d2v(t)dt2=0,(2.12)p1:d2v1(t)dt2+[δ+4εcos⁡(t)2−ϕv0(t)2−2ε]v0(t)=0,(2.13)p2:ε4×103[−102ε+(φ−100δ)t2+2×102−3φ]cos⁡(2t)−tε2×103(φ−102δ)sin⁡(2t)+ε2cos⁡(4t)3.2×102+(2.1875×10−2−2.5×10−2t2)ε2+ε24×105[(−1.8×103φ+6×104δ)t2+1.8×103φ−1.2×105δ]+(φ−102δ/3)(φ−102δ)t48×105.Here,
the initial condition at this boundary can be determined by the boundary
condition (see Section 1).With
the effective initial approximation foru0 from the initial conditions to (2.11), we construct u0:(2.14)v0(t)=0.1.Then,
solving (2.12), (2.13), we have(2.15)v1(t)=(ϕ200−δ2)t2−ε2+εcos⁡(2t)2,v2(t)=14×103ε[−102ε+(−102δ+φ)t2+2×102δ−3φ]cos⁡(2t)−12×103εt(−102δ+φ)sin⁡(2t)+13.2×102ε2cos⁡(4t)+12.4×106(5.25×104−6×104t2)ε2+12.4×106[(−1.8×103φ+6×104δ)t2−1.2×105δ+1.8×103φ]ε+18×105(−102δ+φ)(φ−1023δ)t4.So substituting (2.14)-(2.15) into (2.10)
gives the solution in the following form:(2.16)v(t)=0.1+p[(φ2×102−δ2)t2−ε2+εcos⁡(2t)2]+p2{14×103ε[−102ε+(−102δ+φ)t2+2×102δ−3ψ]cos⁡(2t)−12×103εt(−102δ+φ)sin⁡(2t)+13.2×102ε2cos⁡(4t)+12.4×106(5.25×104−6×104t2)ε2+12.4×106[(−1.8×103φ+6×104δ)t2−1.2×105δ+1.8×103φ]ε+18×105(−102δ+φ)(φ−1023δ)t4}+⋯.As p→1,v(t)→x(t) in (2.16), the solution for the displacement
function of vibrations in a
parametrically excited oscillator with strong cubic negative nonlinearity
is given. From (2.14)-(2.15), it can be realized
that HPM requires only three steps to get more accurate results rather than
other analytical methods.


## 2.1. Homotopy-Perturbation Method (HPM)

### 2.1.1. Background

In this paper, we apply HPM [22–37] for the
solution of the Mathieu-Duffing
equation (1.1). In order to demonstrate how this method works, let us consider
the following nonlinear differential equation:(2.1)A(u)−f(r)=0,r∈Ω,subject to the
boundary  conditions  of(2.2)B(u,∂u∂r)=0,r∈Γ,where
A is a general differential operator, B is a boundary operator, f(r)
is a known analytic functional, and Γ is the boundary of the domain Ω.The
operatorA can generally be divided into two parts L and
N, where L is linear, whereas N is nonlinear. Therefore, (2.1)
can be rewritten as(2.3)L(u)+N(u)−f(r)=0.Homotopy-perturbation
structure is shown to be as the following equation [22–37]:(2.4)H(v,p)=(1−p)[L(v)−L(u0)]+p[L(v)+N(v)−f(r)]=0,where(2.5)v(r,p):Ω×[0,1]→R.For
p = 0 and p = 1, (2.4) reduces to the following equations,
respectively:(2.6)H(v,0)=L(v)−L(u0)=0,H(v,1)=A(v)−f(r)=0,where p∈[0,1] is an embedding parameter and u0 is the first approximation that satisfied the
boundary condition.The
process of changes inp from zero to unity is that of v(r,p) changing form u0 to u(r).
We consider v, as the following [22–37]:(2.7)v=v0+pv1+p2v2+⋯,and
the best approximation for solution is [22–37](2.8)u=lim⁡p→1v=v0+v1+v2+⋯.


### 2.1.2. Implementation of HPM

Now we apply homotopy-perturbation (2.4)
to (1.1). We construct a homotopy in the following form:(2.9)(1−p)[d2x(t)dt2]+p{−[δ+2εcos⁡(2t)]x(t)+ϕx(t)3}=0.Obviously,
(2.9) becomes linear if p = 0, when p = 1, it becomes the original nonlinear
one. So the variation of p from zero to unity makes the equation to change to nonlinear
from linear one.According
to HPM, we assume that the solution of (2.9) can be expressed in a series of p(2.10)v(r)=v0(r)+pv1(r)+p2v2(r)+⋯.Substituting v(r) from (2.10)
into (2.9), after some simplification and substitution and rearranging based on
powers of p-terms, we have(2.11)p0:d2v(t)dt2=0,(2.12)p1:d2v1(t)dt2+[δ+4εcos⁡(t)2−ϕv0(t)2−2ε]v0(t)=0,(2.13)p2:ε4×103[−102ε+(φ−100δ)t2+2×102−3φ]cos⁡(2t)−tε2×103(φ−102δ)sin⁡(2t)+ε2cos⁡(4t)3.2×102+(2.1875×10−2−2.5×10−2t2)ε2+ε24×105[(−1.8×103φ+6×104δ)t2+1.8×103φ−1.2×105δ]+(φ−102δ/3)(φ−102δ)t48×105.Here,
the initial condition at this boundary can be determined by the boundary
condition (see Section 1).With
the effective initial approximation foru0 from the initial conditions to (2.11), we construct u0:(2.14)v0(t)=0.1.Then,
solving (2.12), (2.13), we have(2.15)v1(t)=(ϕ200−δ2)t2−ε2+εcos⁡(2t)2,v2(t)=14×103ε[−102ε+(−102δ+φ)t2+2×102δ−3φ]cos⁡(2t)−12×103εt(−102δ+φ)sin⁡(2t)+13.2×102ε2cos⁡(4t)+12.4×106(5.25×104−6×104t2)ε2+12.4×106[(−1.8×103φ+6×104δ)t2−1.2×105δ+1.8×103φ]ε+18×105(−102δ+φ)(φ−1023δ)t4.So substituting (2.14)-(2.15) into (2.10)
gives the solution in the following form:(2.16)v(t)=0.1+p[(φ2×102−δ2)t2−ε2+εcos⁡(2t)2]+p2{14×103ε[−102ε+(−102δ+φ)t2+2×102δ−3ψ]cos⁡(2t)−12×103εt(−102δ+φ)sin⁡(2t)+13.2×102ε2cos⁡(4t)+12.4×106(5.25×104−6×104t2)ε2+12.4×106[(−1.8×103φ+6×104δ)t2−1.2×105δ+1.8×103φ]ε+18×105(−102δ+φ)(φ−1023δ)t4}+⋯.As p→1,v(t)→x(t) in (2.16), the solution for the displacement
function of vibrations in a
parametrically excited oscillator with strong cubic negative nonlinearity
is given. From (2.14)-(2.15), it can be realized
that HPM requires only three steps to get more accurate results rather than
other analytical methods.


## 2.1.1. Background

In this paper, we apply HPM [22–37] for the
solution of the Mathieu-Duffing
equation (1.1). In order to demonstrate how this method works, let us consider
the following nonlinear differential equation:(2.1)A(u)−f(r)=0,r∈Ω,subject to the
boundary  conditions  of(2.2)B(u,∂u∂r)=0,r∈Γ,where
A is a general differential operator, B is a boundary operator, f(r)
is a known analytic functional, and Γ is the boundary of the domain Ω.The
operatorA can generally be divided into two parts L and
N, where L is linear, whereas N is nonlinear. Therefore, (2.1)
can be rewritten as(2.3)L(u)+N(u)−f(r)=0.Homotopy-perturbation
structure is shown to be as the following equation [22–37]:(2.4)H(v,p)=(1−p)[L(v)−L(u0)]+p[L(v)+N(v)−f(r)]=0,where(2.5)v(r,p):Ω×[0,1]→R.For
p = 0 and p = 1, (2.4) reduces to the following equations,
respectively:(2.6)H(v,0)=L(v)−L(u0)=0,H(v,1)=A(v)−f(r)=0,where p∈[0,1] is an embedding parameter and u0 is the first approximation that satisfied the
boundary condition.The
process of changes inp from zero to unity is that of v(r,p) changing form u0 to u(r).
We consider v, as the following [22–37]:(2.7)v=v0+pv1+p2v2+⋯,and
the best approximation for solution is [22–37](2.8)u=lim⁡p→1v=v0+v1+v2+⋯.


## 2.1.2. Implementation of HPM

Now we apply homotopy-perturbation (2.4)
to (1.1). We construct a homotopy in the following form:(2.9)(1−p)[d2x(t)dt2]+p{−[δ+2εcos⁡(2t)]x(t)+ϕx(t)3}=0.Obviously,
(2.9) becomes linear if p = 0, when p = 1, it becomes the original nonlinear
one. So the variation of p from zero to unity makes the equation to change to nonlinear
from linear one.According
to HPM, we assume that the solution of (2.9) can be expressed in a series of p(2.10)v(r)=v0(r)+pv1(r)+p2v2(r)+⋯.Substituting v(r) from (2.10)
into (2.9), after some simplification and substitution and rearranging based on
powers of p-terms, we have(2.11)p0:d2v(t)dt2=0,(2.12)p1:d2v1(t)dt2+[δ+4εcos⁡(t)2−ϕv0(t)2−2ε]v0(t)=0,(2.13)p2:ε4×103[−102ε+(φ−100δ)t2+2×102−3φ]cos⁡(2t)−tε2×103(φ−102δ)sin⁡(2t)+ε2cos⁡(4t)3.2×102+(2.1875×10−2−2.5×10−2t2)ε2+ε24×105[(−1.8×103φ+6×104δ)t2+1.8×103φ−1.2×105δ]+(φ−102δ/3)(φ−102δ)t48×105.Here,
the initial condition at this boundary can be determined by the boundary
condition (see Section 1).With
the effective initial approximation foru0 from the initial conditions to (2.11), we construct u0:(2.14)v0(t)=0.1.Then,
solving (2.12), (2.13), we have(2.15)v1(t)=(ϕ200−δ2)t2−ε2+εcos⁡(2t)2,v2(t)=14×103ε[−102ε+(−102δ+φ)t2+2×102δ−3φ]cos⁡(2t)−12×103εt(−102δ+φ)sin⁡(2t)+13.2×102ε2cos⁡(4t)+12.4×106(5.25×104−6×104t2)ε2+12.4×106[(−1.8×103φ+6×104δ)t2−1.2×105δ+1.8×103φ]ε+18×105(−102δ+φ)(φ−1023δ)t4.So substituting (2.14)-(2.15) into (2.10)
gives the solution in the following form:(2.16)v(t)=0.1+p[(φ2×102−δ2)t2−ε2+εcos⁡(2t)2]+p2{14×103ε[−102ε+(−102δ+φ)t2+2×102δ−3ψ]cos⁡(2t)−12×103εt(−102δ+φ)sin⁡(2t)+13.2×102ε2cos⁡(4t)+12.4×106(5.25×104−6×104t2)ε2+12.4×106[(−1.8×103φ+6×104δ)t2−1.2×105δ+1.8×103φ]ε+18×105(−102δ+φ)(φ−1023δ)t4}+⋯.As p→1,v(t)→x(t) in (2.16), the solution for the displacement
function of vibrations in a
parametrically excited oscillator with strong cubic negative nonlinearity
is given. From (2.14)-(2.15), it can be realized
that HPM requires only three steps to get more accurate results rather than
other analytical methods.


## 3. Numerical Solution

In order to
provide a basis for the purpose of comparison of analysis of parametrically
excited vibrations of an oscillator with strong cubic negative nonlinearity
with no explicit solution, in this section, the solution of the governing
equation of the system (1.1) by using RK iterative method is presented. The
method can also be
used for other second-order differential equations.


### 3.1. Runge-Kutta's Method (RK)

As shown by (1.1), the solution of vibration
of parametrically excited oscillator with strong cubic negative nonlinearity is
to solve a second-order differential equation with a set of certain boundary
conditions. For this case, (1.1) can be written in the following general form:(3.1)x¨=f(t,x,x˙)=−(δ+2εcos⁡2t)x+ϕx.For such a boundary value problem given by boundary
condition defined in Section 4, some numerical methods have been developed [38, 39]. Here we apply the fourth-order RK algorithm to solve (3.1) subject to the
given boundary conditions. RK iterative formulae for the second-order
differential equation are [39](3.2)x˙i+1=x˙i+Δt6(k1+2k2+2k3+k4),xi+1=xi+Δt[x˙i+Δt6(k1+k2+k3)],where Δt is the increment of the time
and k1,k2,k3, and k4 are determined from the following formulae:(3.3)k1=f(ti,xi,x˙i),k2=f(ti+Δt2,xi+Δt2x˙i,x˙i+Δt2k1),k3=f(ti+Δt2,xi+Δt2x˙i+14Δt2k1,x˙i+Δt2k2),k4=f(ti+Δt,xi+Δtx˙i+12Δt2k2,x˙i+Δtk3).The
numerical solution starts from the boundary at the initial time, where the first
value of the displacement function and its first-order derivative is determined
from initial condition (see Section 1). Then, with a small time increment (Δt), the displacement function and its first-order
derivative at the new position can be obtained using (3.2). This process
continues to the end of time.


## 3.1. Runge-Kutta's Method (RK)

As shown by (1.1), the solution of vibration
of parametrically excited oscillator with strong cubic negative nonlinearity is
to solve a second-order differential equation with a set of certain boundary
conditions. For this case, (1.1) can be written in the following general form:(3.1)x¨=f(t,x,x˙)=−(δ+2εcos⁡2t)x+ϕx.For such a boundary value problem given by boundary
condition defined in Section 4, some numerical methods have been developed [38, 39]. Here we apply the fourth-order RK algorithm to solve (3.1) subject to the
given boundary conditions. RK iterative formulae for the second-order
differential equation are [39](3.2)x˙i+1=x˙i+Δt6(k1+2k2+2k3+k4),xi+1=xi+Δt[x˙i+Δt6(k1+k2+k3)],where Δt is the increment of the time
and k1,k2,k3, and k4 are determined from the following formulae:(3.3)k1=f(ti,xi,x˙i),k2=f(ti+Δt2,xi+Δt2x˙i,x˙i+Δt2k1),k3=f(ti+Δt2,xi+Δt2x˙i+14Δt2k1,x˙i+Δt2k2),k4=f(ti+Δt,xi+Δtx˙i+12Δt2k2,x˙i+Δtk3).The
numerical solution starts from the boundary at the initial time, where the first
value of the displacement function and its first-order derivative is determined
from initial condition (see Section 1). Then, with a small time increment (Δt), the displacement function and its first-order
derivative at the new position can be obtained using (3.2). This process
continues to the end of time.


## 4. Results and Discussion

In this
study, the Mathieu-Duffing equation has
been solved by using HPM and RK method. The
results shown in Figures1–6 indicate that
the HPM experiences a high accuracy. In addition, in comparison with RK method and
other analytical methods, a considerable reduction of the volume of the
calculation can be seen in HPM. It can be approved that HPM is powerful in finding
analytical solutions for a wide class of nonlinear problems.The result of HPM and RK atφ = 2, ε = 0.01, x(0)=0.1, x˙(0)=0.
(a) Time history diagram of x. (b) Time history diagram of x˙.
(a)(b)Figure 2
The result of HPM and RK,x˙(t)  versus  x(t) at φ = 2, ε = 0.01, x(0)=0.1, x˙(0)=0.Figure 3
Accuracy of the results shown in Figure1.Figure 4
Accuracy of the results shown in Figure2.Figure 5
Accuracy of Figure1 versus ε (0.01<ε<1).Figure 6
Accuracy of Figure1 versus φ (2<φ<4).Figures1(a) and 1(b) illustrate the time history diagram of the displacement and
velocity, respectively. Figure 2 illustrates the velocity versus the displacement
obtained from the proposed method. Figure 1 as well as Figure 2 are obtained
for φ = 2, ε = 0.01. These figures show obviously the excellent agreement between HPM and RK method. Figure 2 also shows the residual loop; the area of loops
corresponds to the loss of energy in each cycle. It illustrates that the increase of
cycle numbers reduces the loss of energy.For further verification, the accuracy of Figure1 is shown in Figure 3 which
shows that increasing the time gives rise to reduce the accuracy of displacement as well as the
velocity. Figure 4 confirms the accuracy of results depicted in Figure 2. Figure 5
shows the effects of the variation of the parameter ε at t = 7 and φ = 2 on the accuracy of HPM. From this figure, it is clear that the
change of the ε from 0.01 to 1 reduces the accuracy up to 10%. So, it can be seen that the best approximation
for high accurate results can be obtained in the range of ε less than 0.1. From Figure 6, it is evident that at t = 7 and
ε = 0.01, the variation of φ from 2 to 4 reduces the accuracy up to 95%.


## 5. Conclusion

In this survey,
the HPM has been employed to analyze the parametrically
excited vibrations of an oscillator with strong cubic negative nonlinearity. The results obtained from
this method have been compared with those obtained from numerical method using
RK algorithm. This comparison shows excellent agreement between the two
methods.The presented
scheme provides concise and straightforward solution to approach reliable
results, and it overcomes the difficulties that have been arisen in
conventional methods. Also, HPM does not require small parameters, so the
limitation of the conventional perturbation method could be eliminated.Solution  of
the  Mathieu-Duffing equation shows that accuracy of the
results is affected by the variation of the parametersε and φ considerably. In other words, increasing of
the ε and φ results in a reduction in the accuracy.



---

*Source: 956170-2008-09-08.xml*