- Email: [email protected]
PERIODIC SOLUTIONS OF AN UNSYMMETRIC OSCILLATOR INCLUDING A COMPREHENSIVE STUDY OF THEIR STABILITY CHARACTERISTICS
Journal of Sound and Vibration (1996) 192(5), 959–976
PERIODIC SOLUTIONS OF AN UNSYMMETRIC OSCILLATOR INCLUDING A COMPREHENSIVE STUDY OF THEIR STABILITY CHARACTERISTICS P. D L. N. V Department of Mechanical Engineering and Material Science, Duke University, Durham, NC 27708, U.S.A.
J. J. W U.S. Army Research Office, Research Triangle Park, NC 27709, U.S.A. (Received 13 June 1994, and in final form 15 September 1995) The periodic solutions of a non-linear oscillator with unsymmetric restoring force and harmonic excitation are studied by means of harmonic balance using an arbitrary number of modes in the assumed solution. Comparisons between the approximate solution for two modes are made with respect to both one mode and numerical solutions. Different stability criteria have been used in a comparative analysis which shows that higher order criteria not only give more accurate results but also respect the correct sequence of flips and folds without a significant increase in computational cost. 7 1996 Academic Press Limited
1. INTRODUCTION
The loss of stability of cyclic behavior, as a parameter is slowly varied, is a universal problem in physical sciences. For low order systems described by ordinary differential equations, the generic routes to instability are relatively well defined, and limited in number. In the context of non-linear oscillators, an archetypal unsymmetric oscillator has been the subject of recent investigations: x¨ + bx˙ + x − x 2 = F sin vt,
(1)
where numerical simulation has shown extremely subtle bifurcation behavior, especially associated with the transition to unbounded motion [1, 2]. A number of approximate analytical studies have also been conducted [3, 4, 5]. This paper extends previous analytical work based on the method of harmonic balance. Two main issues are (i) how accurate the solutions are and (ii) how efficient their computation is. Furthermore, the determination of the stability characteristics of periodic solutions is very important. Although these issues are generally assessed in terms of quantitative accuracy (using numerical integration as a basis), the stability characteristics can also be judged by certain qualitative features. Specifically, Thompson [1] has shown that approximate analytical methods may violate the sequence of bifurcations (flip and folds) encountered as a system parameter is slowly changed. In this paper, a comprehensive stability survey is conducted in an attempt to resolve some of the issues outlined above. Periodic solutions to equation (1) are obtained using an arbitrary number of harmonics computed from the general non-linear equations for the 959 0022–460X/96/200959 + 18 $18.00/0
7 1996 Academic Press Limited
. .
960
Fourier coefficients using sophisticated numerical techniques. The local stability of these solutions is then determined using a variety of techniques, which are then compared and judged according to both quantitative and qualitative criteria. 2. COMPUTING PERIODIC SOLUTIONS USING n-MODE HARMONIC BALANCE
Although the focus is on the unsymmetric oscillator (1), the harmonic balance [6–8] procedure is exemplified for the general ordinary differential equation: x˙ = f(x, t),
(2)
where f: DURp × R+:Rp is a periodic function with period T, assuming for the periodic solution the truncated Fourier series n
x˜ = a0 + s (ak cos kvt + bk sin kvt).
(3)
k=1
The periodic function F(t) = f(x˜ (t), t) can be expanded [9] as n
F(t) = A0 + s (Ak cos kvt + Bk sin kvt)
(4)
k=1
and a system of non-linear equations for unknowns a0 , ak , bk , k = 1, n is obtained: A0 = 0,
Ak − kvbk = 0,
Bk + kvak = 0,
k = 1, n .
(5)
Applying the above procedure to the non-linear forced oscillator (1), the number of ordinary differential equations (2) is p = 2 and the following system of algebraic equations is obtained: n
a0 − a02 − 12 s (ak2 + bk2 ) = 0, k=1
(1 − 2a0 − k 2v 2 )ak + bkvbk − 12
s
−bkvak + (1 − 2a0 − k 2v 2 )bk − 12
s
(ai aj − bi bj ) − 12
i,j e 1,i + j = k
s i,j e 1,i + j = k
(ai aj + bi bj ) = 0,
i,j e 1,=i − j= = k
(ai bj + aj bi ) −
s
(aj bi − ai bj ) = Fd1k ,
i,j e i − j = k
k = 1, n , (6) where a0 , ak and bk , k = 1, n , are the Fourier coefficients of the assumed solution. The one-mode approximation (n = 1) has been studied by Virgin [3]. In this paper a two-mode approximation (n = 2) is studied extensively, starting from the following system of algebraic equations: a0 − a02 − 12 (a12 + b12 ) − 12 (a22 + b22 ) = 0, (1 − 2a0 − v 2 )a1 + bvb1 − a1 a2 − b1 b2 = 0, −bva1 + (1 − 2a0 − v 2 )b1 − a1 b2 + a2 b1 = F, (1 − 2a0 − 4v 2 )a2 + 2bvb2 − 12 (a12 − b12 ) = 0, −2bva2 + (1 − 2a0 − 4v 2 )b2 − a1 b1 = 0,
(7)
961
which is numerically solved using a Newton–Raphson method with appropriate path-following capability (based on continuation [10]; see Appendix C). The complete set of real roots (maximum six) has been obtained and their correct number has been found using complex algebra, resulting in a 12th order polynomial equivalent to the system of equations (7) (for which the authors are grateful to Ben Noble). Throughout the paper the damping is fixed at 0·1. A typical time series and phase projection are shown in Figures 1(a) and (b), respectively. These responses are generated for F = 0·1 and v = 0·5 and hence have significant superharmonic component, and the accuracy of the two-mode analysis is apparent. The next few figures show how the response depends on the two main parameters; the forcing magnitude F and frequency v. The maximum peak-to-peak amplitude is presented as a function of v in Figure 2. The maximum displacement of the time series in Figure 1(a) can be located appropriately in Figures 2(c) and (d). A further two mode results are shown in Figure 3. A spurious solution representing oscillation about the unstable equilibrium point has also been found. An alternative and convenient response measure which is analogous to a Poincare´ section is obtained by setting t = 0 in equations (3). For the two-mode analysis, xp = a 0 + a 1 + a 2 ,
x˙p = v(b1 + 2b2 ).
(8)
Solving the relevant set of algebraic equations and placing the coefficients into the above equations gives the results of Figures 4 and 5. A comparison between the one-mode and two-mode results for v = 0·85 is shown in Figure 4, while a summary of two mode results for various forcing frequencies is presented in Figure 5.
Figure 1. A typical steady state oscillation comparing one-mode (- - -), two-mode (· · ·) and numerical integration (----); v = 0·5, F = 0·1: (a) time series; (b) phase projection.
962
. .
Figure 2. Amplitude–frequency response diagrams comparing the one-mode and two-mode solutions: (a) one mode, F = 0·05; (b) two modes, F = 0·05; (c) one mode, F = 0·1; (d) two modes, F = 0·1.
3. LOCAL STABILITY OF PERIODIC SOLUTIONS
The results of the previous section include both stable and unstable solutions. It is important to determine which of these solutions are stable and hence physically realizable. Local stability is studied comparatively through different criteria either analytically or numerically and the results are compared with fully numerical results reported by Thompson [1]. Steady state period-one oscillations typically become unstable via ‘‘flip’’ and ‘‘fold’’ bifurcations. The flip, also termed a period-doubling bifurcation, signals the appearance of an n = 2 subharmonic oscillation, which is often the start of a sequence leading to chaos and occurs when one of the system characteristic multipliers penetrates the unit circle in the complex plane at −1 on the real axis [11]. The fold bifurcation, also termed a saddle node, is often associated with the amplitude resonant jump and hysteresis and is defined by a characteristic multiplier penetrating the unit circle at +1 on the real axis. Of the many methods to analyze stability (e.g., Routh–Hurwitz), Floquet theory provides a good starting point. Considering a small perturbation j around the periodic solution x0 of the second order differential equation x¨ = f(x, x˙ , t; l),
(9)
the following variational equation is obtained: j = rj + bj ,
(10)
963
Figure 3. Amplitude–frequency response diagrams using the two-mode analysis: (a) F = 0·135; (b) an enlarged view; (c) F = 0·15; (d) F = 0·2.
Figure 4. A Poincare´ section comparing the one mode and two mode results for v = 0·85: (a) displacement versus forcing amplitude; (b) velocity versus forcing amplitude; - - -, one-mode solution; ----, two-mode solution.
. .
964
Figure 5. Poincare´ sections based on the two-mode analysis: (a) and (c), v = 0·05; (b) and (d), v = 0·1; (e) and (g), v = 0·4; (f) and (h), v = 0·6.
where r=
1f (x , x˙ , t; l), 1x 0 0
b=
From Floquet theory, j = emth,
1f (x , x˙ , t; l). 1x˙ 0 0
965
where h is a periodic function with the same period, the following linear ordinary differential equation is obtained: h¨ = (b − 2m)h˙ + (r + mb − m 2 )h. Following Hayashi [6], the characteristic exponent m can be real or imaginary but not complex and h has period T or 2T. Using the truncated expansions m
m
b = a0 + s (ak cos kvt + bk sin kvt),
r = r0 + s (rk cos kvt + sk sin kvt), k=1
k=1
n
h = a0 + s l=1
0
al cos
1
lvt lvt + bl sin , 2 2
the following linear system is obtained: m
n
0 = (r0 + ma0 − m 2 )a0 + 12 s s k=1 l=1
$0
rk+mak −
zcv
1 0
1%
lv lv b a + sk + mbk + ak bl , 2 k l 2
2k − l = 0
−
01 jv 2
2
aj = (a0 − 2m)
m
jv b + (r0 + ma0 − m 2 )aj + ejk (rk + mak )a0 2 j
n
+ 12 s s k=1 l=1
zcv
$0
lv lv b a + −sk − mbk + ak bl 2 k l 2
$0
lv lv b a + sk + mbk + ak bl , 2 k l 2
rk + mak +
1 0
1%
2k + l = j m
n
+ 12 s s k=1 l=1
zcv
rk + mak −
1 0
1%
=2k − l= = j
−
01 jv 2
2
bj = −(a0 − 2m)
m
n
+ 12 s s k=1 l=1
zcv
jv a + (r0 + ma0 − m 2 )bj + ejk (sk + mbk )a0 2 j
$0
sk + mbk −
1 0
1%
lv lv a a + rk + mak + bk bl 2 k l 2
2k + l = j m
n
+ 12 s s sgn (2k − l) k=1 l=1
zcv
$0
sk + mbk +
1 0
1%
lv lv a a + −rk − mak + bk bl , 2 k l 2
=2k − l= = j
j = 1, n ,
ejk =
6
1, 0,
1 E 2k = j E m, otherwise.
(11)
. .
966
In the case of the unsymmetric oscillator: m
n
0 = −(1 − 2a0 + mb + m 2 )a0 + s s [ak al + bk bl ], k=1 l=1
zcv 2k = l = 0
01
−
jv 2
2
aj = −(b + 2m)
m
jv b − (1 − 2a0 + mb + m 2 )aj + 2ejk ak a0 2 j
n
m
n
+ s s [ak al − bk bl ] + s s [ak al + bk bl ], k=1 l=1
k=1 l=1
zcv
zcv
= 2k − l = = j
2k + l = j
01
−
jv 2
2
bj = (b + 2m)
m
jv a − (1 − 2a0 + mb + m 2 )bj + 2ejk bk a0 2 j
n
m
n
+ s s [bk al + ak bl ] + s s sgn (2k − l)[bk al − ak bl ], k=1 l=1
k=1 l=1
zcv
zcv =2k − l= = j
2k + l = j
j = 1, n ,
(12)
where ak and bk are the coefficients in the periodic solution (3). The determinant of the above systems should be zero such that they have a non-trivial solution resulting in truncation of Hill’s infinite determinant as the stability criterion [12]. Although both flips and folds can be computed from the above criteria, the folds have been determined using the vertical tangency rule as has been presented by Virgin [3]. In this paper, higher order criteria (i.e., n = 3 and n = 4) have been tested and the result compared with the well-known first order (n = 1) and improved (n = 2) criteria presented in Hayashi [6]. 4. DISCUSSION
The stability boundary for the flip bifurcation (in the v − F space) is defined for both the first order and the improved analysis by solving the system of equations (7) with the corresponding criterion attached (see Appendix B). In Figure 6 are shown loci of flip bifurcations plotted in the v − F plane for the first order stability analysis. Figure 6(a) is based on the one-mode assumption and Figure 6(b) is based on the two-mode assumption. The improved stability criteria are shown in Figure 7(b), together with a comparison with the first order stability result in Figure 7(a). Both of these figures are based on a two-mode assumed solution. The flip bifurcations are also shown in Figure 8, but now projected into the Poincare´ section versus forcing frequency plane. These results are based on the two-mode, improved stability analysis. The fold bifurcation is defined by 1F 1a 1r 1r = 0 \ 0 = 2a \ 1 = 2a \ 2 = 2a. 1a0 1F 1F 1F
(13)
967
Figure 6. Loci of flip bifurcations in the F–v plane: (a) one mode, first order stability analysis; (b) two modes, first order stability analysis.
Figure 7. Loci of flip bifurcations in the F–v plane: (a) two modes, first order stability analysis; (b) two modes, second order stability analysis. - - -, equation (14); ----, equation (15).
Figure 8. Loci of the flip bifurcations plotted against the response (Poincare´ displacement) using two modes and the second order stability analysis: (a)–(c) equation (B2); (d) equation (B3).
. .
968
Figure 9. Loci of the fold bifurcations in the v–F plane: (b) is a detail of (a).
The above conditions are equivalent to the condition that the determinant of the Jacobian of the system (7) vanishes. Starting with system (7), one can use the last two equations to eliminate a2 and b2 . Then the following system of three equations is obtained: a0 − a02 −
r12 r14 − = 0, 2 2 8(A2 + B22 )
A1 b1 − B1 a1 −
r12 A2 b1 + B2 a1 = F, 2 A22 + B22
A1 a1 + B1 b1 −
r12 A2 a1 − B2 b1 = 0. 2 A22 + B22
(14)
By squaring the second and the third equations and adding them, a new equation in z and r1 only is obtained. Therefore the Jacobian of the following reduced equivalent system has to vanish:
a0 − a02 − (A12 + B12 )r12 + r14
r12 r14 1 − = 0, 2 8 A22 + B22
B1 B2 − A1 A2 r16 1 + = F 2, A22 + B22 4 A22 + B22
(15)
where A1 = 1 − 2a0 − v 2,
B1 = bv,
A2 = 1 − 2a0 − 4v 2,
B2 = 2bv.
The fold bifurcations are shown in Figure 9, based on the two-mode solution. There are now two cusps which relate to vertical tangencies in the primary and superharmonic peaks. T 1 Fold bifurcation; F = 0·1 Method
Frequency
Amplitude
1 2, 3, 7 4, 5, 6, 8
0·795 0·771 0·796
0·711 0·712 0·765
969
Higher order criteria (n = 3 and n = 4) have also been analyzed without a significant increase in the computation time due to the matrix bandedness of the linear system of equations (11) and (12). Their motivation lies in the fact that using the same number of modes in the perturbation as in the periodic solution actually ignores the equations and the terms in which the first modes in the perturbation interacts with the last modes in the periodic solution. In order to compare the first order and improved stability criteria with the higher order ones (containing n e 2 modes in the perturbation), a robust engineering criterion is used: Given a certain level of forcing, at what frequency does the steady state solution first lose its stability, what is the type of instability mechanism and what is the amplitude of the response when it occurs? Amplitude (A) versus frequency (v) curves are traced and the stability criterion is recorded. The first point at which instability occurs is recorded for increasing v (typically a fold) and, respectively, decreasing v (typically a flip) in Tables 1–6. The tables contain results from several stability approaches for different steady state solution approximations: (1) numerical simulation (both the solution and local stability are determined numerically); (2) one-mode steady state solution—one mode in the periodic part of the perturbation (first order stability, Hayashi); (3) one-mode steady state solution—two mode in the periodic part of the perturbation (improved criterion, Hayashi); (4) two-mode steady state solution—two-mode in the periodic part of the perturbation (improved criterion, Hayashi); (5) two-mode steady state solution—three-mode in the periodic part of the perturbation (higher order analysis); (6) two-mode steady state solution; four-mode in the periodic part of the perturbation (higher order analysis); (7) numerical integration of the variational equation based on one-mode steady state solution; (8) numerical integration of the variational equation based on two-mode steady state solution. Having established the stability boundaries, and from the comparative analysis of different stability criteria, it is now possible to return to the steady state responses and identify the stable branches. In Figure 10 is shown a typical set of solutions using the Poincare´ sectioning and plotting as a function of forcing amplitude F. Starting from the rest state (i.e., F = xp = x˙p = 0), the response grows until the first flip is encountered prior to the first fold. Since the response curves and stability boundaries are multi-valued and exist in a higher order space, it is difficult to relate intersections between the response and stability boundary simply by following one of the parameters. However, once the system parameters are set the curves are followed simultaneously, thus enabling the appropriate encountering of flips and folds as the solution evolves. Also shown in Figure 10 are results from direct numerical simulation. In Figure 11 is shown a comparison between one-mode, two-mode and numerical simulation in the v–A plane. The two-mode solution is very close to the simulation except for the appearance of a small n = 3 superharmonic peak in the vicinity of v = 0·3. It is also worth mentioning that for this specific F = 0·135, a number of subtle responses are observed and hence the vertical tangencies associated with superharmonic resonant peak are shown to occur for the analytical results and not for the simulation results. In effect, they are very close to the cusp point at v 2 0·45 and F = 0·135 in Figure 9. In general, numerical–analytical correlation is very good. The numerical results are plotted separately (in Figure 11(c)), since they are almost indistinguishable from the two mode analysis when superimposed. Apart from the quantitative evaluation, the stability analysis can be assessed qualitatively by the correct sequence of bifurcations (flip and folds). Thus, the characteristic multipliers can move on a circle inside the unit circle or on the real axis such
. .
970
T 2 Flip bifurcation; F = 0·1 Method
Frequency
Amplitude
1 2 3 4 5 6 7 8
0·852 0·946 — 0·869 0·840 0·839 — —
1·223 0·982 — 1·179 1·221 1·225 — —
T 3 Fold bifurcation; F = 0·135 Method
Frequency
Amplitude
1 2, 3, 7 4, 5, 6, 8
0·741 0·704 0·747
0·778 0·779 0·673
T 4 Flip bifurcation; F = 0·135 Method
Frequency
Amplitude
1 2 3 4 5 6 7 8
0·922 0·994 — 0·928 0·916 0·915 — 0·859
1·159 0·968 — 1·138 1·160 1·162 — 1·242
T 5 Fold bifurcation; F = 0·2 Method
Frequency
Amplitude
1 2, 3, 7 4, 5, 6, 8
0·402 0·577 0·349
0·744 0·861 0·678
T 6 Flip bifurcation; F = 0·2 Method
Frequency
Amplitude
1 2 3 4 5 6 7 8
1·028 1·082 — 1·026 1·020 1·019 — 0·998
1·070 0·929 — 1·068 1·076 1·079 — 1·119
971
Figure 10. A summary of steady state responses indicating stability for v = 0·6 using Poincare´ sections. W, flip bifurcation; w, fold bifurcation. (a) F versus xp ; (b) F versus x˙p ; — — —, numerical.
that their product is less than one as the system is dissipative [1]. Also, if one characteristic multiplier leaves the unit circle through −1 (flip bifurcation), the other characteristic multiplier approaches the origin but cannot pass through it, and if the oscillator has lost its stability through a period-doubling bifurcation it can restabilize only through a reverse period-doubling. The same observations are valid for the fold bifurcation, showing that
Figure 11. A summary of the steady state response indicating stability for F = 0·135, representing peak-to-peak amplitude versus frequency: (a) one-mode solution; (b) two-mode solution; (c) numerical simulation (stable solution).
972
. .
Figure 12. Stable and unstable solutions using Poincare´ sections for v = 0·4: (a) one-mode solution/one-mode perturbation; (b) one-mode solution/two-mode perturbation; (c) two-mode solution/two-mode perturbation; (d) two-mode solution/three-mode perturbation; (e) two-mode solution/four-mode perturbation; (f) zoom of (e).
flip and fold bifurcations should occur only in pairs as the periodic solution is path-followed. The sequencing of flips and folds has been studied for various frequencies (0·4, 0·5, 0·6, 0·75, 0·85, 0·95) showing the same pattern, but only v = 0·4 and v = 0·6 are presented in Figures 12 and 13. The one-mode solution and first order stability analysis (a) and the two-mode solution and improved stability analysis (c) do not satisfy the correct sequence of flips and folds, but an extra mode in the perturbation (i.e., one-mode solution/two-mode perturbation or two-mode solution/three-mode perturbation) respects the correct sequence. A higher order stability analysis does not improve the stability and in one case even the sequence of flips/folds is not observed, which confirms Hamdan and Burton’s remarks on Duffing’s equation [13].
973
Figure 13. Stable and unstable solutions using Poincare´ sections for v = 0·6: (a) one-mode solution/one-mode perturbation; (b) one-mode solution/two-mode perturbation; (c) two-mode solution/two-mode perturbation; (d) two-mode solution/three-mode perturbation; (e) two-mode solution/four-mode perturbation.
5. CONCLUSIONS
Accurate solutions of the period one steady state solutions can be computed by including an arbitrary number of terms in equations (3) and (6). A second order harmonic was included, as an example, since the non-linearity is quadratic and earlier results based on a single harmonic were not satisfactory in the low frequency range when compared with numerical simulation. This suggested inclusion of a term to capture superharmonic effects. Furthermore, the focus of many practical applications in dynamics is behavior associated with primary resonance, and in softening spring oscillators such as equation (1) the resonant peak tends to bend toward lower frequencies. The harmonic balance method has been used to determine approximate solutions. One mode periodic solution gives good accuracy for v q 1 and reasonable results for 0·8 Q v Q 1 (Figures 10 and 11) [3, 1]. By using a two-mode expansion the results can be used for a larger range of frequencies, i.e., 0·3 Q v. Accuracy is improved even for frequencies 0·7 Q v Q 1 because the second mode influences the first mode. A second peak in the frequency response appears at about v = 0·5 (Figure 3). One can predict that including a third harmonic will capture the peak that appears at about v = 0·33. The use of the improved stability condition has a strong effect on the stability analysis of the solutions. For example, there exists a small restabilizing domain which has not been
. .
974
predicted by the use of the first order criterion (see Poincare´ sections for v = 0·6, in Figure 10, for large F). The restabilizing portion has been predicted numerically by Thompson [1]. Fold bifurcations can be followed even to zero frequency with the exception of one curve which could be traced only to 0·08; an ill-conditioning problem is encountered in numerical simulations as well [1, 3, 4]. Also, comparisons between different stability criteria show that a higher order criterion containing an extra mode in the perturbation gives the most accurate results without a significant increase in the computation cost and a correct sequence of flips and folds. ACKNOWLEDGMENT
This work has been supported by grant number DAAL03-89-D-0002 from the U.S. Army Research Office. REFERENCES 1. J. M. T. T 1989 Proceedings of the Royal Society London A 421, 195–225. Chaotic phenomena triggering the escape from a potential well. 2. S. F and J. M. T. T 1991 Computer Methods in Applied Mechanics and Engineering 89, 381–394. Geometrical concepts and computational techniques of nonlinear dynamics. 3. L. N. V 1988 Journal of Sound and Vibration 126, 157–165. On the harmonic response of an oscillator with unsymmetric restoring force. 4. E. R, A. R-L and M. G. V 1992 in From Phase Transitions to Chaos, (G. Gyorgyi, I. Kondor, L. Savari and T. Tel, editors). London: World Scientific. Oscillations and onset of chaos in a driven nonlinear oscillator with possible escape to infinity. 5. J. L. S 1995 Nonlinear Dynamics 7, 11–35. Variable coefficient harmonic balance for periodically forced nonlinear oscillators. 6. C. H 1985 Nonlinear Oscillations in Physical Systems. Princeton, New Jersey: Princeton University Press. 7. D. W. J and P. S 1987 Nonlinear Ordinary Differential Equations. Oxford: Clarendon Press. 8. W. S-S and J. R 1993 Physica D 66, 368–380. Bifurcation phenomena in a nonlinear oscillator: approximate analytical studies versus computer simulation results. 9. F. H. L and X. X. W 1987 International Journal of Non-Linear Mechanics 22(2), 89–98. Fast Galerkin method and its application to determine periodic solutions of non-linear oscillators. 10. R. S 1991 International Journal of Bifurcation and Chaos 1(1), 3–11. Tutorial on continuation. 11. L. N. V and K. D. M 1994 Journal of Sound and Vibration 169, 699–703. On the behavior of characteristic multipliers through a period-doubling sequence. 12. S. T. N and G. R. H 1982 Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics 49, 217–223. A generalized Hill’s method for the stability analysis of parametrically excited dynamic systems. 13. M. N. H and T. D. B 1993 Journal of Sound and Vibration 166, 255–266. On the steady state response and stability of non-linear oscillators using harmonic balance. 14. B. N 1964 Numerical Methods. Oliver & Boyd.
APPENDIX A: DERIVATION OF THE NON-LINEAR ALGEBRAIC SYSTEM
As an example one of the coefficients of the expansion (4) is computed: Ak =
2 T
g
T
0
F(t) cos kvt dt =
2 T
g
T
0
f(x˜ (t), t) cos kvt dt.
(A1)
975
For the oscillator under consideration:
f(x(t), t) =
$
%
x2 , −(x1 − x12 ) − bx2 + F sin vt
x=
$ %
%$x1 x = . x2 x˙
(A2)
The coefficients in the solution expansion (3) are
ak =
$ %
a1,k , k = 0, n , a2,k
bk =
$ %
b1,k , b2,k
k = 1, n .
(A3)
Introducing equations (A3), (A2) and (A1) in equation (5) and integrating one obtains a2,k − kvb1,k = 0,
−ba2,k − a1,k + 2a1,0 a1,k + 12
s
(a1,i a1,j − b1,i b1,j )
i,j e 1,i + j = k
− 12
s
(a1,i a1,j + b1,i b1,j ) − kvb2,k = 0.
i,j e 1,=i − j= = k
From the computation of Bk , b2,k + kva1,k = 0. Eliminating a2,k and b2,k and omitting the subscript 1 in a1,k and b1,k , the second set of equations (6) is obtained. The other equations are obtained in a similar manner.
APPENDIX B: STABILITY CRITERIA
The first order stability criterion is given by [6]
01 0 1
0 1 01
2 G n v G u0 + d 2 − − unc u1s − 2 n d G 2 2 G G G q 0, 2 G n G v 2 u0 + d − u1s + 2 n d + unc G 2 G 2
where unc = −an , uns = −bn , n = 1, 2.
(B1)
. .
976
The improved stability condition, i.e., a second order approximation, is [6] 2 G v v u1s − d u1c − u2c G u0 + d2 − 2 − u1c 2 G 2 G v v u1s + d u0 + d 2 − + u1c u1s + u2s G 2 2 G G 3v u1c − u2c u1s + u2s u0 + d 2 − G 2 G G −u1s + u2s u1c + u2c 3vd G G
01
01
0 1
G G G G u1c + u2c Gq 0 G G −3vd G G 2G 3v u0 + d 2 − G 2 G −u1s + u2s
2
0 1
(B2) for the ‘‘first’’ unstable region, and
n
u0 + d 2 2u1s 2u1c
n
u1s u1c u0 + d 2 − v 2 − u2c u2s − 2vd q0 u0 + d 2 − v 2 + u2c u2s + 2vd
(B3)
for the ‘‘second’’ unstable region. APPENDIX C: NUMERICAL PATH-FOLLOWING
The harmonic balance method and the stability analysis involve solving systems of non-linear equations of the type f(x; l) = 0,
(C1)
where x $ Rn is the unknown and l $ R is the parameter. Continuation methods based on the arc-length method are used for path following the solution. The predictor is based on the tangent for starting the path following and after the first step uses the secant [10, 2]. The corrector step involves solving a system of non-linear equations consisting of the initial equation and a constraint z(Dx)T(Dx) + (Dl)2 = Ds,
(C2)
where Dx is the step in the unknown x, Dl is the step in the parameter l, and Ds is the arc-length parameter. The algebraic system with x and l as unknowns is solved using the well-known Newton’s method [14] with backtracking as line search.
PERIODIC SOLUTIONS OF AN UNSYMMETRIC OSCILLATOR INCLUDING A COMPREHENSIVE STUDY OF THEIR STABILITY CHARACTERISTICS P. D L. N. V Department of Mechanical Engineering and Material Science, Duke University, Durham, NC 27708, U.S.A.
J. J. W U.S. Army Research Office, Research Triangle Park, NC 27709, U.S.A. (Received 13 June 1994, and in final form 15 September 1995) The periodic solutions of a non-linear oscillator with unsymmetric restoring force and harmonic excitation are studied by means of harmonic balance using an arbitrary number of modes in the assumed solution. Comparisons between the approximate solution for two modes are made with respect to both one mode and numerical solutions. Different stability criteria have been used in a comparative analysis which shows that higher order criteria not only give more accurate results but also respect the correct sequence of flips and folds without a significant increase in computational cost. 7 1996 Academic Press Limited
1. INTRODUCTION
The loss of stability of cyclic behavior, as a parameter is slowly varied, is a universal problem in physical sciences. For low order systems described by ordinary differential equations, the generic routes to instability are relatively well defined, and limited in number. In the context of non-linear oscillators, an archetypal unsymmetric oscillator has been the subject of recent investigations: x¨ + bx˙ + x − x 2 = F sin vt,
(1)
where numerical simulation has shown extremely subtle bifurcation behavior, especially associated with the transition to unbounded motion [1, 2]. A number of approximate analytical studies have also been conducted [3, 4, 5]. This paper extends previous analytical work based on the method of harmonic balance. Two main issues are (i) how accurate the solutions are and (ii) how efficient their computation is. Furthermore, the determination of the stability characteristics of periodic solutions is very important. Although these issues are generally assessed in terms of quantitative accuracy (using numerical integration as a basis), the stability characteristics can also be judged by certain qualitative features. Specifically, Thompson [1] has shown that approximate analytical methods may violate the sequence of bifurcations (flip and folds) encountered as a system parameter is slowly changed. In this paper, a comprehensive stability survey is conducted in an attempt to resolve some of the issues outlined above. Periodic solutions to equation (1) are obtained using an arbitrary number of harmonics computed from the general non-linear equations for the 959 0022–460X/96/200959 + 18 $18.00/0
7 1996 Academic Press Limited
. .
960
Fourier coefficients using sophisticated numerical techniques. The local stability of these solutions is then determined using a variety of techniques, which are then compared and judged according to both quantitative and qualitative criteria. 2. COMPUTING PERIODIC SOLUTIONS USING n-MODE HARMONIC BALANCE
Although the focus is on the unsymmetric oscillator (1), the harmonic balance [6–8] procedure is exemplified for the general ordinary differential equation: x˙ = f(x, t),
(2)
where f: DURp × R+:Rp is a periodic function with period T, assuming for the periodic solution the truncated Fourier series n
x˜ = a0 + s (ak cos kvt + bk sin kvt).
(3)
k=1
The periodic function F(t) = f(x˜ (t), t) can be expanded [9] as n
F(t) = A0 + s (Ak cos kvt + Bk sin kvt)
(4)
k=1
and a system of non-linear equations for unknowns a0 , ak , bk , k = 1, n is obtained: A0 = 0,
Ak − kvbk = 0,
Bk + kvak = 0,
k = 1, n .
(5)
Applying the above procedure to the non-linear forced oscillator (1), the number of ordinary differential equations (2) is p = 2 and the following system of algebraic equations is obtained: n
a0 − a02 − 12 s (ak2 + bk2 ) = 0, k=1
(1 − 2a0 − k 2v 2 )ak + bkvbk − 12
s
−bkvak + (1 − 2a0 − k 2v 2 )bk − 12
s
(ai aj − bi bj ) − 12
i,j e 1,i + j = k
s i,j e 1,i + j = k
(ai aj + bi bj ) = 0,
i,j e 1,=i − j= = k
(ai bj + aj bi ) −
s
(aj bi − ai bj ) = Fd1k ,
i,j e i − j = k
k = 1, n , (6) where a0 , ak and bk , k = 1, n , are the Fourier coefficients of the assumed solution. The one-mode approximation (n = 1) has been studied by Virgin [3]. In this paper a two-mode approximation (n = 2) is studied extensively, starting from the following system of algebraic equations: a0 − a02 − 12 (a12 + b12 ) − 12 (a22 + b22 ) = 0, (1 − 2a0 − v 2 )a1 + bvb1 − a1 a2 − b1 b2 = 0, −bva1 + (1 − 2a0 − v 2 )b1 − a1 b2 + a2 b1 = F, (1 − 2a0 − 4v 2 )a2 + 2bvb2 − 12 (a12 − b12 ) = 0, −2bva2 + (1 − 2a0 − 4v 2 )b2 − a1 b1 = 0,
(7)
961
which is numerically solved using a Newton–Raphson method with appropriate path-following capability (based on continuation [10]; see Appendix C). The complete set of real roots (maximum six) has been obtained and their correct number has been found using complex algebra, resulting in a 12th order polynomial equivalent to the system of equations (7) (for which the authors are grateful to Ben Noble). Throughout the paper the damping is fixed at 0·1. A typical time series and phase projection are shown in Figures 1(a) and (b), respectively. These responses are generated for F = 0·1 and v = 0·5 and hence have significant superharmonic component, and the accuracy of the two-mode analysis is apparent. The next few figures show how the response depends on the two main parameters; the forcing magnitude F and frequency v. The maximum peak-to-peak amplitude is presented as a function of v in Figure 2. The maximum displacement of the time series in Figure 1(a) can be located appropriately in Figures 2(c) and (d). A further two mode results are shown in Figure 3. A spurious solution representing oscillation about the unstable equilibrium point has also been found. An alternative and convenient response measure which is analogous to a Poincare´ section is obtained by setting t = 0 in equations (3). For the two-mode analysis, xp = a 0 + a 1 + a 2 ,
x˙p = v(b1 + 2b2 ).
(8)
Solving the relevant set of algebraic equations and placing the coefficients into the above equations gives the results of Figures 4 and 5. A comparison between the one-mode and two-mode results for v = 0·85 is shown in Figure 4, while a summary of two mode results for various forcing frequencies is presented in Figure 5.
Figure 1. A typical steady state oscillation comparing one-mode (- - -), two-mode (· · ·) and numerical integration (----); v = 0·5, F = 0·1: (a) time series; (b) phase projection.
962
. .
Figure 2. Amplitude–frequency response diagrams comparing the one-mode and two-mode solutions: (a) one mode, F = 0·05; (b) two modes, F = 0·05; (c) one mode, F = 0·1; (d) two modes, F = 0·1.
3. LOCAL STABILITY OF PERIODIC SOLUTIONS
The results of the previous section include both stable and unstable solutions. It is important to determine which of these solutions are stable and hence physically realizable. Local stability is studied comparatively through different criteria either analytically or numerically and the results are compared with fully numerical results reported by Thompson [1]. Steady state period-one oscillations typically become unstable via ‘‘flip’’ and ‘‘fold’’ bifurcations. The flip, also termed a period-doubling bifurcation, signals the appearance of an n = 2 subharmonic oscillation, which is often the start of a sequence leading to chaos and occurs when one of the system characteristic multipliers penetrates the unit circle in the complex plane at −1 on the real axis [11]. The fold bifurcation, also termed a saddle node, is often associated with the amplitude resonant jump and hysteresis and is defined by a characteristic multiplier penetrating the unit circle at +1 on the real axis. Of the many methods to analyze stability (e.g., Routh–Hurwitz), Floquet theory provides a good starting point. Considering a small perturbation j around the periodic solution x0 of the second order differential equation x¨ = f(x, x˙ , t; l),
(9)
the following variational equation is obtained: j = rj + bj ,
(10)
963
Figure 3. Amplitude–frequency response diagrams using the two-mode analysis: (a) F = 0·135; (b) an enlarged view; (c) F = 0·15; (d) F = 0·2.
Figure 4. A Poincare´ section comparing the one mode and two mode results for v = 0·85: (a) displacement versus forcing amplitude; (b) velocity versus forcing amplitude; - - -, one-mode solution; ----, two-mode solution.
. .
964
Figure 5. Poincare´ sections based on the two-mode analysis: (a) and (c), v = 0·05; (b) and (d), v = 0·1; (e) and (g), v = 0·4; (f) and (h), v = 0·6.
where r=
1f (x , x˙ , t; l), 1x 0 0
b=
From Floquet theory, j = emth,
1f (x , x˙ , t; l). 1x˙ 0 0
965
where h is a periodic function with the same period, the following linear ordinary differential equation is obtained: h¨ = (b − 2m)h˙ + (r + mb − m 2 )h. Following Hayashi [6], the characteristic exponent m can be real or imaginary but not complex and h has period T or 2T. Using the truncated expansions m
m
b = a0 + s (ak cos kvt + bk sin kvt),
r = r0 + s (rk cos kvt + sk sin kvt), k=1
k=1
n
h = a0 + s l=1
0
al cos
1
lvt lvt + bl sin , 2 2
the following linear system is obtained: m
n
0 = (r0 + ma0 − m 2 )a0 + 12 s s k=1 l=1
$0
rk+mak −
zcv
1 0
1%
lv lv b a + sk + mbk + ak bl , 2 k l 2
2k − l = 0
−
01 jv 2
2
aj = (a0 − 2m)
m
jv b + (r0 + ma0 − m 2 )aj + ejk (rk + mak )a0 2 j
n
+ 12 s s k=1 l=1
zcv
$0
lv lv b a + −sk − mbk + ak bl 2 k l 2
$0
lv lv b a + sk + mbk + ak bl , 2 k l 2
rk + mak +
1 0
1%
2k + l = j m
n
+ 12 s s k=1 l=1
zcv
rk + mak −
1 0
1%
=2k − l= = j
−
01 jv 2
2
bj = −(a0 − 2m)
m
n
+ 12 s s k=1 l=1
zcv
jv a + (r0 + ma0 − m 2 )bj + ejk (sk + mbk )a0 2 j
$0
sk + mbk −
1 0
1%
lv lv a a + rk + mak + bk bl 2 k l 2
2k + l = j m
n
+ 12 s s sgn (2k − l) k=1 l=1
zcv
$0
sk + mbk +
1 0
1%
lv lv a a + −rk − mak + bk bl , 2 k l 2
=2k − l= = j
j = 1, n ,
ejk =
6
1, 0,
1 E 2k = j E m, otherwise.
(11)
. .
966
In the case of the unsymmetric oscillator: m
n
0 = −(1 − 2a0 + mb + m 2 )a0 + s s [ak al + bk bl ], k=1 l=1
zcv 2k = l = 0
01
−
jv 2
2
aj = −(b + 2m)
m
jv b − (1 − 2a0 + mb + m 2 )aj + 2ejk ak a0 2 j
n
m
n
+ s s [ak al − bk bl ] + s s [ak al + bk bl ], k=1 l=1
k=1 l=1
zcv
zcv
= 2k − l = = j
2k + l = j
01
−
jv 2
2
bj = (b + 2m)
m
jv a − (1 − 2a0 + mb + m 2 )bj + 2ejk bk a0 2 j
n
m
n
+ s s [bk al + ak bl ] + s s sgn (2k − l)[bk al − ak bl ], k=1 l=1
k=1 l=1
zcv
zcv =2k − l= = j
2k + l = j
j = 1, n ,
(12)
where ak and bk are the coefficients in the periodic solution (3). The determinant of the above systems should be zero such that they have a non-trivial solution resulting in truncation of Hill’s infinite determinant as the stability criterion [12]. Although both flips and folds can be computed from the above criteria, the folds have been determined using the vertical tangency rule as has been presented by Virgin [3]. In this paper, higher order criteria (i.e., n = 3 and n = 4) have been tested and the result compared with the well-known first order (n = 1) and improved (n = 2) criteria presented in Hayashi [6]. 4. DISCUSSION
The stability boundary for the flip bifurcation (in the v − F space) is defined for both the first order and the improved analysis by solving the system of equations (7) with the corresponding criterion attached (see Appendix B). In Figure 6 are shown loci of flip bifurcations plotted in the v − F plane for the first order stability analysis. Figure 6(a) is based on the one-mode assumption and Figure 6(b) is based on the two-mode assumption. The improved stability criteria are shown in Figure 7(b), together with a comparison with the first order stability result in Figure 7(a). Both of these figures are based on a two-mode assumed solution. The flip bifurcations are also shown in Figure 8, but now projected into the Poincare´ section versus forcing frequency plane. These results are based on the two-mode, improved stability analysis. The fold bifurcation is defined by 1F 1a 1r 1r = 0 \ 0 = 2a \ 1 = 2a \ 2 = 2a. 1a0 1F 1F 1F
(13)
967
Figure 6. Loci of flip bifurcations in the F–v plane: (a) one mode, first order stability analysis; (b) two modes, first order stability analysis.
Figure 7. Loci of flip bifurcations in the F–v plane: (a) two modes, first order stability analysis; (b) two modes, second order stability analysis. - - -, equation (14); ----, equation (15).
Figure 8. Loci of the flip bifurcations plotted against the response (Poincare´ displacement) using two modes and the second order stability analysis: (a)–(c) equation (B2); (d) equation (B3).
. .
968
Figure 9. Loci of the fold bifurcations in the v–F plane: (b) is a detail of (a).
The above conditions are equivalent to the condition that the determinant of the Jacobian of the system (7) vanishes. Starting with system (7), one can use the last two equations to eliminate a2 and b2 . Then the following system of three equations is obtained: a0 − a02 −
r12 r14 − = 0, 2 2 8(A2 + B22 )
A1 b1 − B1 a1 −
r12 A2 b1 + B2 a1 = F, 2 A22 + B22
A1 a1 + B1 b1 −
r12 A2 a1 − B2 b1 = 0. 2 A22 + B22
(14)
By squaring the second and the third equations and adding them, a new equation in z and r1 only is obtained. Therefore the Jacobian of the following reduced equivalent system has to vanish:
a0 − a02 − (A12 + B12 )r12 + r14
r12 r14 1 − = 0, 2 8 A22 + B22
B1 B2 − A1 A2 r16 1 + = F 2, A22 + B22 4 A22 + B22
(15)
where A1 = 1 − 2a0 − v 2,
B1 = bv,
A2 = 1 − 2a0 − 4v 2,
B2 = 2bv.
The fold bifurcations are shown in Figure 9, based on the two-mode solution. There are now two cusps which relate to vertical tangencies in the primary and superharmonic peaks. T 1 Fold bifurcation; F = 0·1 Method
Frequency
Amplitude
1 2, 3, 7 4, 5, 6, 8
0·795 0·771 0·796
0·711 0·712 0·765
969
Higher order criteria (n = 3 and n = 4) have also been analyzed without a significant increase in the computation time due to the matrix bandedness of the linear system of equations (11) and (12). Their motivation lies in the fact that using the same number of modes in the perturbation as in the periodic solution actually ignores the equations and the terms in which the first modes in the perturbation interacts with the last modes in the periodic solution. In order to compare the first order and improved stability criteria with the higher order ones (containing n e 2 modes in the perturbation), a robust engineering criterion is used: Given a certain level of forcing, at what frequency does the steady state solution first lose its stability, what is the type of instability mechanism and what is the amplitude of the response when it occurs? Amplitude (A) versus frequency (v) curves are traced and the stability criterion is recorded. The first point at which instability occurs is recorded for increasing v (typically a fold) and, respectively, decreasing v (typically a flip) in Tables 1–6. The tables contain results from several stability approaches for different steady state solution approximations: (1) numerical simulation (both the solution and local stability are determined numerically); (2) one-mode steady state solution—one mode in the periodic part of the perturbation (first order stability, Hayashi); (3) one-mode steady state solution—two mode in the periodic part of the perturbation (improved criterion, Hayashi); (4) two-mode steady state solution—two-mode in the periodic part of the perturbation (improved criterion, Hayashi); (5) two-mode steady state solution—three-mode in the periodic part of the perturbation (higher order analysis); (6) two-mode steady state solution; four-mode in the periodic part of the perturbation (higher order analysis); (7) numerical integration of the variational equation based on one-mode steady state solution; (8) numerical integration of the variational equation based on two-mode steady state solution. Having established the stability boundaries, and from the comparative analysis of different stability criteria, it is now possible to return to the steady state responses and identify the stable branches. In Figure 10 is shown a typical set of solutions using the Poincare´ sectioning and plotting as a function of forcing amplitude F. Starting from the rest state (i.e., F = xp = x˙p = 0), the response grows until the first flip is encountered prior to the first fold. Since the response curves and stability boundaries are multi-valued and exist in a higher order space, it is difficult to relate intersections between the response and stability boundary simply by following one of the parameters. However, once the system parameters are set the curves are followed simultaneously, thus enabling the appropriate encountering of flips and folds as the solution evolves. Also shown in Figure 10 are results from direct numerical simulation. In Figure 11 is shown a comparison between one-mode, two-mode and numerical simulation in the v–A plane. The two-mode solution is very close to the simulation except for the appearance of a small n = 3 superharmonic peak in the vicinity of v = 0·3. It is also worth mentioning that for this specific F = 0·135, a number of subtle responses are observed and hence the vertical tangencies associated with superharmonic resonant peak are shown to occur for the analytical results and not for the simulation results. In effect, they are very close to the cusp point at v 2 0·45 and F = 0·135 in Figure 9. In general, numerical–analytical correlation is very good. The numerical results are plotted separately (in Figure 11(c)), since they are almost indistinguishable from the two mode analysis when superimposed. Apart from the quantitative evaluation, the stability analysis can be assessed qualitatively by the correct sequence of bifurcations (flip and folds). Thus, the characteristic multipliers can move on a circle inside the unit circle or on the real axis such
. .
970
T 2 Flip bifurcation; F = 0·1 Method
Frequency
Amplitude
1 2 3 4 5 6 7 8
0·852 0·946 — 0·869 0·840 0·839 — —
1·223 0·982 — 1·179 1·221 1·225 — —
T 3 Fold bifurcation; F = 0·135 Method
Frequency
Amplitude
1 2, 3, 7 4, 5, 6, 8
0·741 0·704 0·747
0·778 0·779 0·673
T 4 Flip bifurcation; F = 0·135 Method
Frequency
Amplitude
1 2 3 4 5 6 7 8
0·922 0·994 — 0·928 0·916 0·915 — 0·859
1·159 0·968 — 1·138 1·160 1·162 — 1·242
T 5 Fold bifurcation; F = 0·2 Method
Frequency
Amplitude
1 2, 3, 7 4, 5, 6, 8
0·402 0·577 0·349
0·744 0·861 0·678
T 6 Flip bifurcation; F = 0·2 Method
Frequency
Amplitude
1 2 3 4 5 6 7 8
1·028 1·082 — 1·026 1·020 1·019 — 0·998
1·070 0·929 — 1·068 1·076 1·079 — 1·119
971
Figure 10. A summary of steady state responses indicating stability for v = 0·6 using Poincare´ sections. W, flip bifurcation; w, fold bifurcation. (a) F versus xp ; (b) F versus x˙p ; — — —, numerical.
that their product is less than one as the system is dissipative [1]. Also, if one characteristic multiplier leaves the unit circle through −1 (flip bifurcation), the other characteristic multiplier approaches the origin but cannot pass through it, and if the oscillator has lost its stability through a period-doubling bifurcation it can restabilize only through a reverse period-doubling. The same observations are valid for the fold bifurcation, showing that
Figure 11. A summary of the steady state response indicating stability for F = 0·135, representing peak-to-peak amplitude versus frequency: (a) one-mode solution; (b) two-mode solution; (c) numerical simulation (stable solution).
972
. .
Figure 12. Stable and unstable solutions using Poincare´ sections for v = 0·4: (a) one-mode solution/one-mode perturbation; (b) one-mode solution/two-mode perturbation; (c) two-mode solution/two-mode perturbation; (d) two-mode solution/three-mode perturbation; (e) two-mode solution/four-mode perturbation; (f) zoom of (e).
flip and fold bifurcations should occur only in pairs as the periodic solution is path-followed. The sequencing of flips and folds has been studied for various frequencies (0·4, 0·5, 0·6, 0·75, 0·85, 0·95) showing the same pattern, but only v = 0·4 and v = 0·6 are presented in Figures 12 and 13. The one-mode solution and first order stability analysis (a) and the two-mode solution and improved stability analysis (c) do not satisfy the correct sequence of flips and folds, but an extra mode in the perturbation (i.e., one-mode solution/two-mode perturbation or two-mode solution/three-mode perturbation) respects the correct sequence. A higher order stability analysis does not improve the stability and in one case even the sequence of flips/folds is not observed, which confirms Hamdan and Burton’s remarks on Duffing’s equation [13].
973
Figure 13. Stable and unstable solutions using Poincare´ sections for v = 0·6: (a) one-mode solution/one-mode perturbation; (b) one-mode solution/two-mode perturbation; (c) two-mode solution/two-mode perturbation; (d) two-mode solution/three-mode perturbation; (e) two-mode solution/four-mode perturbation.
5. CONCLUSIONS
Accurate solutions of the period one steady state solutions can be computed by including an arbitrary number of terms in equations (3) and (6). A second order harmonic was included, as an example, since the non-linearity is quadratic and earlier results based on a single harmonic were not satisfactory in the low frequency range when compared with numerical simulation. This suggested inclusion of a term to capture superharmonic effects. Furthermore, the focus of many practical applications in dynamics is behavior associated with primary resonance, and in softening spring oscillators such as equation (1) the resonant peak tends to bend toward lower frequencies. The harmonic balance method has been used to determine approximate solutions. One mode periodic solution gives good accuracy for v q 1 and reasonable results for 0·8 Q v Q 1 (Figures 10 and 11) [3, 1]. By using a two-mode expansion the results can be used for a larger range of frequencies, i.e., 0·3 Q v. Accuracy is improved even for frequencies 0·7 Q v Q 1 because the second mode influences the first mode. A second peak in the frequency response appears at about v = 0·5 (Figure 3). One can predict that including a third harmonic will capture the peak that appears at about v = 0·33. The use of the improved stability condition has a strong effect on the stability analysis of the solutions. For example, there exists a small restabilizing domain which has not been
. .
974
predicted by the use of the first order criterion (see Poincare´ sections for v = 0·6, in Figure 10, for large F). The restabilizing portion has been predicted numerically by Thompson [1]. Fold bifurcations can be followed even to zero frequency with the exception of one curve which could be traced only to 0·08; an ill-conditioning problem is encountered in numerical simulations as well [1, 3, 4]. Also, comparisons between different stability criteria show that a higher order criterion containing an extra mode in the perturbation gives the most accurate results without a significant increase in the computation cost and a correct sequence of flips and folds. ACKNOWLEDGMENT
This work has been supported by grant number DAAL03-89-D-0002 from the U.S. Army Research Office. REFERENCES 1. J. M. T. T 1989 Proceedings of the Royal Society London A 421, 195–225. Chaotic phenomena triggering the escape from a potential well. 2. S. F and J. M. T. T 1991 Computer Methods in Applied Mechanics and Engineering 89, 381–394. Geometrical concepts and computational techniques of nonlinear dynamics. 3. L. N. V 1988 Journal of Sound and Vibration 126, 157–165. On the harmonic response of an oscillator with unsymmetric restoring force. 4. E. R, A. R-L and M. G. V 1992 in From Phase Transitions to Chaos, (G. Gyorgyi, I. Kondor, L. Savari and T. Tel, editors). London: World Scientific. Oscillations and onset of chaos in a driven nonlinear oscillator with possible escape to infinity. 5. J. L. S 1995 Nonlinear Dynamics 7, 11–35. Variable coefficient harmonic balance for periodically forced nonlinear oscillators. 6. C. H 1985 Nonlinear Oscillations in Physical Systems. Princeton, New Jersey: Princeton University Press. 7. D. W. J and P. S 1987 Nonlinear Ordinary Differential Equations. Oxford: Clarendon Press. 8. W. S-S and J. R 1993 Physica D 66, 368–380. Bifurcation phenomena in a nonlinear oscillator: approximate analytical studies versus computer simulation results. 9. F. H. L and X. X. W 1987 International Journal of Non-Linear Mechanics 22(2), 89–98. Fast Galerkin method and its application to determine periodic solutions of non-linear oscillators. 10. R. S 1991 International Journal of Bifurcation and Chaos 1(1), 3–11. Tutorial on continuation. 11. L. N. V and K. D. M 1994 Journal of Sound and Vibration 169, 699–703. On the behavior of characteristic multipliers through a period-doubling sequence. 12. S. T. N and G. R. H 1982 Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics 49, 217–223. A generalized Hill’s method for the stability analysis of parametrically excited dynamic systems. 13. M. N. H and T. D. B 1993 Journal of Sound and Vibration 166, 255–266. On the steady state response and stability of non-linear oscillators using harmonic balance. 14. B. N 1964 Numerical Methods. Oliver & Boyd.
APPENDIX A: DERIVATION OF THE NON-LINEAR ALGEBRAIC SYSTEM
As an example one of the coefficients of the expansion (4) is computed: Ak =
2 T
g
T
0
F(t) cos kvt dt =
2 T
g
T
0
f(x˜ (t), t) cos kvt dt.
(A1)
975
For the oscillator under consideration:
f(x(t), t) =
$
%
x2 , −(x1 − x12 ) − bx2 + F sin vt
x=
$ %
%$x1 x = . x2 x˙
(A2)
The coefficients in the solution expansion (3) are
ak =
$ %
a1,k , k = 0, n , a2,k
bk =
$ %
b1,k , b2,k
k = 1, n .
(A3)
Introducing equations (A3), (A2) and (A1) in equation (5) and integrating one obtains a2,k − kvb1,k = 0,
−ba2,k − a1,k + 2a1,0 a1,k + 12
s
(a1,i a1,j − b1,i b1,j )
i,j e 1,i + j = k
− 12
s
(a1,i a1,j + b1,i b1,j ) − kvb2,k = 0.
i,j e 1,=i − j= = k
From the computation of Bk , b2,k + kva1,k = 0. Eliminating a2,k and b2,k and omitting the subscript 1 in a1,k and b1,k , the second set of equations (6) is obtained. The other equations are obtained in a similar manner.
APPENDIX B: STABILITY CRITERIA
The first order stability criterion is given by [6]
01 0 1
0 1 01
2 G n v G u0 + d 2 − − unc u1s − 2 n d G 2 2 G G G q 0, 2 G n G v 2 u0 + d − u1s + 2 n d + unc G 2 G 2
where unc = −an , uns = −bn , n = 1, 2.
(B1)
. .
976
The improved stability condition, i.e., a second order approximation, is [6] 2 G v v u1s − d u1c − u2c G u0 + d2 − 2 − u1c 2 G 2 G v v u1s + d u0 + d 2 − + u1c u1s + u2s G 2 2 G G 3v u1c − u2c u1s + u2s u0 + d 2 − G 2 G G −u1s + u2s u1c + u2c 3vd G G
01
01
0 1
G G G G u1c + u2c Gq 0 G G −3vd G G 2G 3v u0 + d 2 − G 2 G −u1s + u2s
2
0 1
(B2) for the ‘‘first’’ unstable region, and
n
u0 + d 2 2u1s 2u1c
n
u1s u1c u0 + d 2 − v 2 − u2c u2s − 2vd q0 u0 + d 2 − v 2 + u2c u2s + 2vd
(B3)
for the ‘‘second’’ unstable region. APPENDIX C: NUMERICAL PATH-FOLLOWING
The harmonic balance method and the stability analysis involve solving systems of non-linear equations of the type f(x; l) = 0,
(C1)
where x $ Rn is the unknown and l $ R is the parameter. Continuation methods based on the arc-length method are used for path following the solution. The predictor is based on the tangent for starting the path following and after the first step uses the secant [10, 2]. The corrector step involves solving a system of non-linear equations consisting of the initial equation and a constraint z(Dx)T(Dx) + (Dl)2 = Ds,
(C2)
where Dx is the step in the unknown x, Dl is the step in the parameter l, and Ds is the arc-length parameter. The algebraic system with x and l as unknowns is solved using the well-known Newton’s method [14] with backtracking as line search.