#### A Dynamic Branch-Switching Method for Parametrically Excited Systems

Shock and Vibration
1070-9622
A dynamic branch-switching method for parametrically excited systems
A.Y.T. Leung 0
T. Ge 0
0 School of Engineering, University of Manchester , Manchester M13 9PL , UK
The branch-switching algorithm in static is applied to steady state dynamic problems. The governing ordinary differential equations are transformed to nonlinear algebraic equations by means of harmonic balance method using multiple frequency components. The frequency components of the (irrational) nonlinearity of oscillator are obtained by Fast Fourier Transform and Toeplitz Jacobian method (FFT/TJM). All singularities, folds, flips, period doubling and period bubbling, are computed accurately in an analytical manner. Coexisting solutions can be predicted without using initial condition search. The consistence of both stability criteria in time and frequency domains is discussed. A highly nonlinear parametrically excited system is given as example. All connected solution paths are predicted.
1. Introduction
A nonlinear dynamical system can be described
accurately by a set of nonlinear ordinary differential
equations in time variable when the spatial variables
are eliminated by means of the finite element method.
It is well known that a nonlinear dynamical system will
have multiple solutions [9]. The number of solutions
is not known before hand and may change as the
system parameters vary. There are qualitative and
quantitative studies of the solutions of a dynamical system
in its parameter space. For steady state numerical
solutions, there are time integration methods [32],
frequency methods [13,30,39] and time-frequency
methods [1]. The time methods are classified by the ways of
replacing the time derivatives by differences: implicit,
explicit and Runge–Kutta methods. The time methods
are simple to implement and can predict many kinds of
responses (superharmonic, subharmonic and chaotic)
reasonably well. However, at singularities, they do not
converge and the trajectories look like chaotic. For
higher order subharmonic solutions, the time step must
be very small so that it takes exceedingly long time to
reach the steady state and it will be very difficult to
distinguish between round off errors and chaos. All
possible combinations of initial conditions must be tried in a
blanket search manner for a set of system parameters to
find all possible steady state solutions. Some sensitive
and/or unstable solutions can be missed easily.
Alternatively, the frequency methods replace the
single time variable by many frequency parameters and
the differential equations are transformed into
algebraic equations, which are much easier to handle than
the former. Newtonian method is very popular in
solving nonlinear equations. The incremental harmonic
balance (IHB) method [13] linearizes the ordinary
differential equations before using the harmonic balance
method. The resulting incremental algebraic equations
are linear and therefore, the original IHB method can
not study bifurcation and chaos [4]. To find
subharmonic responses the initial trial solution must be
obtained by other means, e.g., time method or the
singular perturbation method. Leung and Fung [16] apply
harmonic balance to the ordinary differential
equations first to obtain a set of nonlinear algebraic
equations which are linearized by the Newtonian method
afterwards. Higher order dynamic period doublings
and chaos are predicted accurately and are proved
to be consistent with those of Ueda [38] by analog
simulation. The frequency methods rely on the
possibility of expanding the nonlinearity in terms of the
Fourier series. They can not handle irrational
functions. Lewandowski [27,28] made the method in a
systematic manner.
The first author and his colleagues have studied the
frequency method extensively. Leung and Fung [17]
extends it to nonlinear steady state vibration of frames
by finite element method. They applied the method
to study dynamic snap through [18], spinning
structures [19], elasto-plastic systems [20], two-harmonic
excitation [21], and dynamic sub-structuring [22].
Nonlinear modal analysis was suggested [26].
Fergusson and Leung [3] gave an explicit formula to find
the harmonic balance matrix using matrix
manipulation. Leung and Chui [5] studied coupled Duffing
oscillators and the nonlinear vibration of von-Karman’s
square plate [14]. Leung and Ge [23] discovered that
the Jacobian matrix in the harmonic balance is
actually a Toeplitz matrix which can be represented by a
few vectors. Computational efficiency is greatly
improved. Together with FFT technique, the method is
called the FFT/TJM (Toeplitz Jacobian Matrix). The
FFT/TJM method was applied to discontinuous
oscillators [6], nonlinear multi-modes of Euler beams [24],
double-period excitation [7], and construction of
invariant torus [8].
The time-frequency methods [13] are essentially
frequency methods except that the Fourier expansion of
nonlinearity is replaced by the Fourier transform of
the nonlinear functions which are evaluated in discrete
time. One of the time-frequency methods is called the
alternating frequency time (AFT) method [13], which
has three major drawbacks. (i) The Fourier series is
simply truncated at a prescribed number of terms and
the leading terms are not necessarily balanced. (ii) The
Jacobian is approximately given by difference quotient
which is not accurate enough to locate the
singularity point. (iii) A cut-off fraction is pre-set to determine
the significant terms to be included in the
computation to reduce the number of unknowns. The method is
not conformed to the spirit of the Liapunov–Schmidt
reduction and important nonlinear phenomena
associated with the neglected terms can be missed easily.
The dynamic branch switching method presented is
also a time-frequency method. As the system
parameters vary, our aims are (i) to follow the existing
solution branch; (ii) to find the singular points; and (iii)
to switch to a new branch when the exiting branch
becomes unstable at a singular point. Our method is
different from the AFT method in that (i) the Fourier
series is not simply truncated but rather averaged so
that the retaining terms are balanced in Galerkin’s
sense; (ii) the Jacobian is also averaged analytically
which is possible for piecewise continuous functions;
(iii) branch switching strategy is adopted when the
original path becomes unstable near the bifurcation
points; and (iv) the stability of the solution is
determined by Floquet–Liapunov theory.
Parametrically excited nonlinear systems have been
studied extensively. Watt and Cartmell [41] studied
an externally loaded parametric oscillator by
multipletime-scale. Szabelski and Warminski [36] investigated
a self-excited system with parametric and external
excitation by harmonic balance using two terms
approximation. Kahraman and Blankenship [12] used
harmonic balance to find the interactions between
commensurate parametric and forcing excitation with
clearance. Sanchez and Nayfeh [35] investigated the
global behavior of a biased nonlinear oscilltor under
external and parametric excitations by varying the bias
parameter and solved the resulting equations by
analog computer, FFT and Poincare maps. We use the
FFT/TJM method to study the parametrically excited
system until chaos occurs.
An equation representing a parametrically excited
buckled beam will be taken as an example. Period
doubling and period bubbling (period doubling and
then period halving) curves are determined accurately.
Sudden change from deterministic to chaotic response
within 0.5% change of a system parameter is predicted
without difficulty. The complete response curve for
a parametrically excited system is given for the first
time.
2. Basic formulation
Consider a scalar ODE in the form:
f (x, x_ , x, , t)
F
g(t) = 0,
(
1
)
where f is an analytical function, t is the time, x is
the unknown response, is the fundamental frequency
of the linearized system, F and g are the amplitude
and time function of external force, respectively. Over
dot represents differentiation with respect to time () =
d()=dt = d()=d in which is a non-dimensional
time factor, = t. Solution of Eq. (
1
) can be
expressed in a Fourier series
1
x( ) = a0 + X
2
i=1
where a0, ai and bi are the unknown Fourier coefficient
and L is the order of subharmonics. If one considers the
solution precise to subharmonic order 3, then L = 3.
It will be shown later that chaotic motion via period
doubling occurs when the periodic solution bifurcates
from L = 4 to L = 8 and bifurcates to higher
order subharmonic at an accelerated rate. Set increments
x, , F and expand Eq. (
1
) in Taylor series near
x0, 0, F0 which satisfy Eq. (
1
),
f (x0 + x, x_ 0 + x_ , x0 + x, 0 +
, )
(F0 + F )g( )
= f0 +
x
(F0 + F )g + HOT = 0,
(
3
)
where HOT denotes the higher order terms. The
incremental equation is
in which
M ( ) 02 x + C( ) 0 x_ + K( ) x
= r( ) + F g( )
q( )
,
,
q( ) = 2M ( ) 0 x + C( ) x_
,
and r( ) is the residual given by
where
(
4
)
or, in partitioned matrix form,
{X} = [a0, a1, : : : , aN , b1, : : : , bN ]T,
{ X} = [ a0, a1, : : : , aN , b1, : : : , bN ]T
2 J oo
6 J co
4
J so
J oc
J cc
J cs
J os 3
J sc 7
J ss 5
X +
8 Q0 9
>< Q =>
>: Qcs >;
8 P0 9 8 R0 9
F < P = + < R = = {0},
: P cs ;
: Rcs ;
r( ) = f (x, x_ , x, , )
F0g( ) ) 0:
The solution of the incremental equation satisfies the
condition of harmonic balance, so that
1
a0 + X
2
i=1
x( ) =
is then reached by the Newtonian method when the
residual of the approximate solution approaches to
zero.
The other step in the IHB is the formation of the
Jacobi matrix for Newtonian iteration by the Galerkin
method
(
5
)
(
6
)
j
L
j
L
j
L
where
j 0C( ) sin
cos
d ,
j 0C( ) sin
sin
d ,
K( )
j2 02M ( ) sin
K( )
j2 02M ( ) sin
+ j 0C( ) cos
cos
d ,
+ j 0C( ) cos
sin
d ,
j
L
j
L
j
L
(
8
)
(
9
)
(
10
)
g( ) cos
g( ) sin
q( ) cos
q( ) sin
r( ) cos
r( ) sin
An improvement is now presented using
windowing and sampling techniques to discretize the
continuous integration in Eqs (
11
) and (
12
). The Galerkin
averaging is achieved most efficiently using the Fast
Fourier Transformation (FFT). The results obtained are
exactly the same as those found by Eqs (
11
) and (
12
).
Compared with the method described in the previous
section, a very small amount of analytical work is
required before iteration. The alternative transformations
from the frequency domain to the time domain are
included in order to make use of the advantage of each
domain. A fast transformation pair of real form based
on the standard FFT algorithm [31] is developed in this
scheme to bridge the response between the time
domain and the frequency domain, the basic algorithm
being provided as follows. The order of subharmonic
is assumed to be L.
a) The backward transformation from a0, ai, bi to its
time domain responses xk
,
k = 1, 2, : : : , M:
(
13
)
cM i = ci ,
k = 1, 2, : : : , N ,
The standard inverse FFT algorithm can solve Eq. (
14
).
b) The forward transformation from xk to a0, ai,
and bi
c0 =
1 MX1 xk,
M
k=0
a0 =
ai =
bi =
=
=
The first and second terms in ai and bi are
corresponding to the complex coefficients ci and cM i from
the results of the standard forward FFT. It should be
noted that the discrete Fourier transform implies
periodicity in both the time and frequency domains. The
correct application of the sampling theorem [31] for
M > 2N + 1, where M is a power of 2, must be
satisfied
with
N
rk = R0 + X Rci cos
2
i=1
2 Lki
M
+ Rsi cos
2 Lki
M
M 1
R0 = X rk,
,
,
k=0
Jocjc =
M2 MX1 Akj ,
k=0
MX1 Akj e i 2 MLki
k=0
M 1
+ X Akj e i
2 Lk(m i)
M
MX1 Akj e i 2 MLki
k=0
M 1
X Akj e i
k=0
2 Lk(m i)
M
,
:
Similar discretizations are also available for Jocjs, Jicjs,
and Jisjs.
Therefore, we can find the solution {Q}, { x},
{ x_ } and { x} accordingly by Eqs (
14
) and (
15
). It
was found [23] that the Jacobian matrix thus obtained
is actually a Toeplitz matrix which can be represented
by a few vectors. Efficient computational method can
be derived.
4. Double pivoting arc-length method
Note that in Eq. (
10
), the number of incremental is
more than the number of equations available due to the
presence of and F . One of the increments must
be fixed so that the other increment can be solved from
Eq. (
10
). Therefore, we can rewrite Eq. (
10
) to
G(X,
X, ,
) =
= {R},
(
24
)
X
b
where G is an augmented matrix with order of 2N + 1
(
20
)
(
21
)
(
22
)
(
23
)
(
28
)
(
29
)
We pivot the active increment on the fastest changed
variable [25] in the last augmentation step. Assume
that the active increment in the current iteration is Xk.
The incremental linear Eqs (
10
) can then be solved for
the 2N + 1 unknowns i, where i are scale factors
whose function is similar to directional derivatives,
xli = i xlk ,
i = 1, 2, : : : , k
1, k + 1, : : : , 2N + 2:
(
27
)
The remaining unknown derivative xk= l is solved
from the constraint Eq. (
27
), which can be rewritten as
The subjective Jacobian {Gk} of Eq. (
25
) will be
regular for all turning points since we have chosen the most
fastest changed variable as active at each augmentation
step. The solution of the iteration are expressed as:
the active increment, which is either or F . If the
active increment vanishes, an improved solution with
diminishing residual vector {R} is guaranteed by the
nature of the Newtonian method. The process is
repeated until the magnitude of the residual vector {R} is
acceptably small and a solution (an equilibrium state)
is obtained. This process is called iteration. On the
other hand, if the value of the active increment is
augmented artificially from a known equilibrium state to
find a new solution, the process is called an
augmentation. An alternative application of augmentation and
iteration is an useful strategy in the IHB scheme. The
method mentioned above is successful in predicting
the dynamic response curves in general. However, for
typical nonlinear stability phenomena, this strategy is
not sufficient since the Jacobian is in the neighborhood
of turning points and bifurcation points. Furthermore,
when the system damping is very small, it may also
fail to converge near the sharp peaks of the response
curve. This imperfection is usually improved by
appending an additional equation to the solution scheme.
A well-known example of solving strategy is the
socalled arc-length method as widely discussed in
nonlinear FEM field by Risk [34], Crisfield [2], Wagner
and Wriggers [40]. The basic idea of such methods is
to append an additional equation to the original set of
Eq. (
10
) constraining the active increments by a given
arc-length. Although this method has been
successful in nonlinear static FEM, it needs some
modifications when applies to dynamic problems [29]. A
double pivoting arc-length method is applied here to avoid
the singular Jacobian at turning points and bifurcation
points.
Denote Xn+1 = and enlarge the vector {X} from
order 2N + 1 to order 2N + 2 to include
Eq. (
24
) can be written as
2N + 2 matrix. The arc-length l is then determined by
normalizing the length of the tangent vector
x2 , : : : ,
l
xN+1
l
which provides an additional equation
x1
l
x1
l
,
2
+
+
5. Stability points and secondary branches
It is well known that the physical response curve will
switch to the new bifurcated branch when the stability
of original branch is lost, or one can only get an
unstable solution with its amplitude growing
monotonically [5]. The cascade diagram, which represents a
sequence of period doubling, is the most celebrated
example of stability loss. Although the algorithm of IHB
can take practically all the harmonic terms into
consideration, it is shown by numerous computation
experiments that subharmonic responses can never be
obtained automatically without using branch-switching
algorithm. The emphasis in the following discussion
will focus on the computational aspects. Our objective
is to develop a numerical scheme which is consistent
with the context of our IHB scheme mentioned above
for bifurcation problems of period doubling and period
halving.
i =
i
k k
II =
X i 0
I
Xi 0 +
i
Xk ,
i = 1, 2, : : : , 2N + 1, i 6= k,
giving the tangent vector along the secondary branch
in which {Xi}0 is the nearest bifurcation points and
superscripts I, II represent the solution branches. The
pivoting of the control increment depends on the
character of bifurcation points. For instance, in case of
symmetry breaking the most effective pivoting could be on
the first element in vector { Xi} while for period
double bifurcation, pivoting on the one of the coefficients
of the associated subharmonic terms will successfully
lead to the solution on the secondary branch. The value
of { Xk} can be approximated as
{ Xk} =
k l,
X k
=
(
30
)
(
31
)
(
32
)
(
33
)
{u_ } =
A(t) {u},
2
[A(t)] = 4
C(t)
M (t)
1
K (t) 3
M (t) 5 :
0
1 Z ti
#
ti 1
period response into consideration. The calculation of
the double period branch switching requires additional
consideration. Near a stability point, the dimension of
equations is enlarged so that subharmonic terms are
accommodated. And then the eigenvalues of the
Jacobian are used to find the umber of possible branches
near the bifurcation point. The number of branches is
determined by the number of zero eigenvalues i of
Eq. (
30
).
The eigenvectors i indicate the directions of the
secondary branches associated with i = 0. The
branch switching is performed by substituting the
factors i in the Eq. (
28
) with the scaled eigenvector
This value is critical for a successful branch switching.
6. Floquet–Liapunov theory
When a solution is obtained, its stability
characteristic can be checked by the Floquet–Liapunov theory.
Then the stability of Eq. (
34
) is checked by
evaluating the eigenvalues of the transition matrix [Q] at the
An approximate method based on this theory was
developed by Hsu [11], and Hsu and Cheng [10] who
discretized the state transition matrix by a series of step
function and uses a single pass numerical integration
to a small perturbation. The basic idea is presented as
follows. When the solution is perturbed by x, the
autonomous linear incremental equation is
M ( ) 02 x + C( ) 0 x_ + K ( ) x = 0,
(
34
)
where
,
,
in which M (t), C(t), K (t) are time dependent
functions. Rewrite Eq. (
34
) as
where {u} = [ x, x]T , and [A(t)] periodic with
period T
Discretize the time domain and replace [A(t)] with step
function:
(
35
)
(
36
)
(
37
)
(
38
)
where the sampling period # = T =M and M is the
number of discrete points in the defined period. The
perturbed system can be approximated by
{u_ i} =
Bi(t) {ui}:
Let [Q] be the matrix solution satisfying
[Q] = [BM ][Q],
Q(0) = [I ]:
end of the specified period from time t = M T to time
t = (M + 1)T . The two eigenvalues are so called
Floquet multipliers. If the absolute magnitudes of the
Floquet multipliers are less than unity, the solution is
stable. For any small perturbation x, its oscillation is
diminishing. The explicit form of [Q] can be written as
For dissipative systems, the solutions are stable when
the moduli of the complex eigenvalues are within the
unit circle. When the system parameter changes,
instability occurs if one of its eigenvalues leaving the unit
circle. That a positive real eigenvalue leaves the unit
circle while the other remains inside predicts a limit
point in folded solution curves. The solution on the
folded part is shown by a shift of the equilibrium
position. When a negative real eigenvalue leaves the unit
circle through 1, the point is known as flip
bifurcation point and an unstable one period solution will
bifurcate into a stable double period solution. This
phenomenon is most commonly met at bifurcation points
and happened in the system we considered below.
When a pair of complex conjugate eigenvalue leaves
the unit circle, it is called Hopf bifurcation which has
been discussed by many authors.
7. Numerical example
As an example, we consider a class of problem
whose excitation appears as a coefficient in the
equation of motion. There is a number of earlier studies on
chaos in parametrically excited systems related to the
forced motion of buckled beams [33,37]. The resulting
equation can be written in terms of x as Fig. 1,
!02x + !0cx_ + [1 + q cos
] sin x
= q cos
:
!02 x + !0c x_ + [1 + q cos ] cos x0 x + r
which is periodic with period T = 2 L. The
number of harmonic terms and residual tolerance used in
each iteration depends on the required order of the
subharmonic solutions. The total number of terms must
be increased when large unbalanced residual terms are
detected. We observe later that the sequence of
period double bifurcation occurs at closely
neighbouring points. The evaluation of bifurcation and branch
switching is possible only if the unbalanced residual
can be kept small. in particular, for period 4 orbits,
we have used 2NMAX + 1 = 57 terms and for period
8 orbit 2NMAX + 1 = 113 terms, which are able to
keep the residual less than 10 6 and 10 8, receptively.
8. Results and discussion
For small values of q (q < 0:28637), the solutions
are attracted only by the stable fixed point x = 0,
x = 0 with two Floquet multipliers less than unit in
period 2 . It should be noted that our method eliminates
the interference of the initial conditions and connect all
possible fixed point of the equation with one pseudo
arc-length curve. So no additional efforts are required
to predict the three coexisting solutions of Eq. (42),
say at q = 0.6. Figure 2 shows the folded curve, there
are three solutions, two stables and one unstable, in
some area. At point a, the original period 1(
2
)
solution loses stability with one of its Floquet multipliers
leaving the unit circle through 1. The stable period
2(
4
) solutions could be obtained by using dynamic
branch switching. The period 1 solution regains its
stability at point b where the period 2 solutions merge
with period 1. A period bubbling is completed. At point
c we have a positive real Floquet multiplier exceeding
unit which relates to a fold or a saddle node
bifurcation. The solution is then attracted by the unstable fixed
point x = , x = 0, tracing along the upper
unstable part of the curve with period 2 . The curve folds
back again at point d where the period 1 solution is
attracted the fixed point x = , x = 0 and all Floquet
multipliers stay in the unit cycle. Period doubling
cascade happens when q reaches 0.2869544. We combine
the arc-length and branch-switching algorithm to trace
the stable branches. Figure 3 shows the period 2
stable branch bifurcating into period 4 . At point f the
branches bifurcates again to 8 , and at point g the
stable solution period bifurcates into 16 (Fig. 4). Further
period double bifurcations on the stable branches are
out of trace due to limited digital precision offered by
the computer.
If the determinant of the Jacobian in Eq. (
10
) takes
the double period into consideration, it is called the
enlarged determinant. The zero points of the enlarged
determinant are correspondent to the points where the
Floquet multipliers pass through the unit cycle. This
consistent relation in the determination of solution
stability is very clear from theories and also verified
by our numerous computational practices. The benefit
from the enlarged determinant check is that additional
calculation to obtain Floquet multipliers of solution is
avoided. So based on this we can use the inverse
iteration to find the position of local bifurcation point. We
Fig. 7c. Spectrum q = 0:287360.
use the Floquet–Liapunov theory at the beginning of
bifurcation to check the stability of the present branch
and then we use the enlarged determinant to determine
the exact local bifurcation point and to switch to new
stable branch. The period doubling cascade is
successfully constructed within a very narrow range of
parameter change (from q = 0:2864 to 0.2873685). The
phase diagram, time histories and spectral at specific
values of q are shown in Figs 5–7, respectively.
9. Conclusion
A numerical method for determining the steady state
forced vibration of a highly non-linear system has
been presented. The method in which an alternative
FFT/TJM algorithm is used, is comparatively easy to
formulate the harmonic balance equations as compared
to that of IHB. An effective algorithm to search for the
specified stable branches and characterise the stability
of solutions has been described. The newly developed
method was applied to a forced parametric excitation
system. A number of well-known nonlinear
phenomena such as fold and period double bifurcation is
obtained. The results generated here are coincided with
that of direct numerical integration. The computer time
is greatly reduced and the response diagrams are neatly
constructed.
Journal of
Engineering
Volume 2014
Volume 2014
Journal of
Sensors
Volume 2014
International Journal of
Hindawi Publishing Corporation
ht p:/ www.hindawi.com
Active
and
Advances in
Civil Engineering
Hindawi Publishing Corporation
ht p:/ www.hindawi.com
Journal of
Robotics
Hindawi Publishing Corporation
ht p:/ www.hindawi.com
Advances
ctronics
Submit
your manuscr ipts
VLSI
Design
Hindawi Publishing Corporation
ht p:/ www.hindawi.com
Hindawi Publishing Corporation
Hindawi Publishing Corporation
Navigation and
Observation
Hindawi Publishing Corporation
ht p:/ www.hindawi.com
Modelling
Sim
ulation
&
Engineering
International Journal of
Engineering
Electrical
and
Computer
International Journal of
Aerospace
Engineering
[1] T.M. Cameron and J.H. Griffin , An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems , ASME, J. of Applied Mechanics 56 ( 1989 ), 149 - 154 .
[2] M.A. Crisfield , A fast incremental/iterative solution procedure that handles 'Snap-through' , Computers and Structures 13 ( 1981 ), 55 - 62 .
[3] N.J. Fergusson and A.Y.T. Leung , Harmonic balance calculations by using matrices , J. of Sound and Vibration 182 ( 1995 ), 559 - 563 .
[4] A.A. Ferri , On the equivalence of the incremental harmonic balance method and the harmonic balance - Newton Raphson method , ASME, J. of Applied Mechanics 53 ( 1986 ), 455 - 457 .
[5] P. Friedmann , C.E. Hammond and T.H. Woo , Efficient numerical treatment of periodic systems with application to stability problem , Int. J. for Numerical Methods in Engineering 11 ( 1977 ), 1117 - 1136 .
[6] T. Ge and A.Y.T. Leung , A Toeplitz Jacobian matrix/FFT method for steady state analysis of discontinuous oscillators , Shock and Vibration 2 ( 1995 ), 205 - 218 .
[7] T. Ge and A.Y.T. Leung , Toeplitz Jacobian method for nonlinear double-periodic excitations , Shock and Vibration 4 ( 1997 ), 351 - 359 .
[8] T. Ge and A.Y.T. Leung , Construction of invariant torus using Toeplitz Jacobian Matrices/Fast Fourier Transform approach , Nonlinear Dynamics 15 ( 1998 ), 283 - 305 .
[9] P.J. Holmes and F.C Moon , Strange attractors and chaos in nonlinear mechanics , ASME, J. of Applied Mechanics 50 ( 1983 ), 1021 - 1032 .
[10] C.S. Hsu and W.H. Cheng, Application of the theory of impulsive parametric excitation and new treatment of general parametric excitation problems , ASME, J. of Applied Mechanics 40 ( 1973 ), 8 - 86 .
[11] C.S. Hsu , Impulsive parametric excitation: theory, ASME , J. of Applied Mechanics 39 ( 1972 ), 551 - 558 .
[12] A. Kahraman and G.W. Blankenship , Iteration between commensurate parametric and forcing excitations in a system with clearance , J. of Sound and Vibration 194 ( 1996 ), 317 - 336 .
[13] S.L. Lau , Y.K. Cheung and S.Y. Wu , A variable parameter incremental method for dynamic instability of linear and nonlinear elastic systems , ASME, J. of Applied Mechanics 49 ( 1982 ), 849 - 853 .
[14] A.Y.T. Leung and S.K. Chui , On the non-linear vibration of von-Karman's square plate by the IHB method , J. of Sound and Vibration 204 ( 1997 ), 239 - 247 .
[15] A.Y.T. Leung and S.K. Chui , Nonlinear vibration of coupled Duffing oscillators by an improved incremental harmonic balance method , J. of Sound and Vibration 181 ( 1995 ), 619 - 633 .
[16] A.Y.T. Leung and T.C. Fung , Construction of chaotic regions , J. of Sound and Vibration 131 ( 3 ) ( 1989 ), 445 - 455 .
[17] A.Y.T. Leung and T.C. Fung , Nonlinear steady state vibration of frames by finite element method , Int. J. for Numerical Methods in Engineering 28 ( 1989 ), 1599 - 1618 .
[18] A.Y.T. Leung and T.C. Fung , Nonlinear steady state vibration and dynamic snap through of arches , J. of Earthquake Engineering and Structural Dynamics 19 ( 1990 ), 409 - 430 .
[19] A.Y.T. Leung and T.C. Fung , Geometrically nonlinear vibration of spinning structures , J. of Sound and Vibration 139 ( 1990 ), 43 - 62 .
[20] A.Y.T. Leung and T.C. Fung , Analytical solutions of elastoplastic systems , J. of Sound and Vibration 142 ( 1990 ), 175 - 182 .
[21] A.Y.T. Leung and T.C. Fung , Bifurcations of an oscillator to two-harmonic excitation , Communications in Applied Numerical Methods 6 ( 1990 ), 573 - 582 .
[22] A.Y.T. Leung and T.C. Fung , Linear-nonlinear dynamic substructures , Int. J. for Numerical Methods in Engineering 31 ( 1991 ), 967 - 985 .
[23] A.Y.T. Leung and T. Ge , Toeplitz Jacobian matrix for nonlinear periodic vibration , ASME, J. of Applied Mechanics 117 ( 1995 ), 709 - 717 .
[24] A.Y.T. Leung and T. Ge , Normal multi-modes of nonlinear Euler beams , J. of Sound and Vibration 202 ( 1997 ), 145 - 160 .
[25] A.Y.T. Leung , Nonlinear natural vibration analysis by selective coefficient increment , Computational Mechanics 5 ( 1989 ), 73 - 80 .
[26] A.Y.T. Leung , Nonlinear modal analysis of frames by the incremental harmonic balance method, Dynamics and Stability of Systems 7 ( 1992 ), 43 - 59 .
[27] R. Lewandowski , Computational formulation for periodic vibration of geometrically nonlinear structures - Part 1: Theoretical background , Int. J. of Solids and Structures 34 ( 1997 ), 1925 - 1947 .
[28] R. Lewandowski , Computational formulation for periodic vibration of geometrically nonlinear structures - Part 2: Numerical strategy and examples , Int. J. of Solids and Structures 34 ( 1997 ), 1949 - 1964 .
[29] M. Marek and I. Shreiber , Chaotic Behavior of Derterministic Dissipative Systems , Cambridge University Press, 1991 .
[30] A.H. Nayfeh and D.T. Mook , Nonlinear Oscillations, Wiley, New York, 1979 .
[31] H.J. Nussbanner , Fast Fourier Transform and Convolution Algorithm , Berlin, Springer, 1982 .
[32] T.S. Parker and O.C. Leon , Practical Numerical Algorithms for Chaotic Systems , Springer Verlag, New York, 1989 .
[33] R.H. Plaut and J.C. Hsieh , Chaos in a mechanism with time delays under parametric and external excitation , J. of Sound and Vibration 114 ( 1987 ), 73 - 90 .
[34] E. Risk , Bifurcation and stability - a numerical approach, Innovative Methods for Nonlinear Problems , W.K. Lin , T. Belytschko and K.C. Park, eds, Pineridge, Swansea, 1984 , pp. 313 - 344 .
[35] N.E. Sanchez and A.H. Nayfeh , Global behaviour of a biased nonlinear oscillator under external and parametric excitation , J. of Sound and Vibration 207 ( 1996 ), 137 - 149 .
[36] K. Szabelski and J. Warminski , Self-excited system vibrations with parametric and external excitations , J. of Sound and Vibration 187 ( 1995 ), 595 - 607 .
[37] W. Szemplinska-Stupnica , R.H. Plaut and J.C. Hsieh , Period doubling and chaos in unsymmetrical structures under parametric excitation , ASME, J. of Applied Mechanics 56 ( 1989 ), 947 - 952 .
[38] Y. Ueda , Steady motions exhibited by Duffing's equation: A picture book of regular and chaotic motions , in: New Approaches to Nonlinear Problems in Dynamics, P.J. Holmes, ed., SIAM , Philadelphia, 1980 , pp. 311 - 322 .
[39] M. Urabe , Galerkin's procedure for nonlinear periodic systems , Archives for Rational Mechanics and Analysis 20 ( 1966 ), 120 - 152 .
[40] W. Wagner and P. Wriggers , A simple method for the calculation of postcritical branches , Engineering Computation 5 ( 1988 ), 103 - 109 .
[41] D. Watt and M. Cartmell , An externally loaded parametric oscillation , J. of Sound and Vibration 170 ( 1994 ), 339 - 364 .
and Volume 2014