- Email: [email protected]
Criteria for chaos of a three-well potential oscillator with homoclinic and heteroclinic orbits
Journal
of Sound and
Vibration (1990) 136(l), 17-34
CRITERIA POTENTIAL
FOR
CHAOS
OSCILLATOR
AND
OF A THREE-WELL WITH
HETEROCLINIC G.
HOMOCLINIC ORBITS
X. Lr
Department of Theoretical and Applied Mechanics AND
F. C.
MOON
Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, New York 14853, U.S.A. (Received 13 July 1988, and in revised form 10 April 1989)
Criteria for chaotic vibrations of a non-linear elastic beam having three stabie and two unstable equilibrium positions are compared. The responses of the system exhibit chaotic behavior in certain regions in the driving frequency and amplitude plane. With the use of the Melnikov technique, necessary conditions for chaos are derived based on homoclinic and heteroclinic bifurcations of the system. Numerical simulations of the stable and unstable manifolds, Poincare maps, Lyapunov exponents and fractal basin boundaries are performed. The results indicate that both homoclinic and heteroclinic bifurcations are necessary conditions for chaos; thus, at least two necessary conditions must be satisfied for any chaotic responses. However, the final steady state of motion for many initial conditions becomes unpredictable (in the sense of fractal basin boundaries) if at least one homoclinic or heteroclinic bifurcation occurs. Finally, an experimental criterion for chaotic motion, and an heuristic criterion based on classical methods in non-linear dynamics are also obtained and compared with the bifurcation criteria, as a function of forcing amplitude and frequency.
1. INTRODUCTION
non-linear oscillator is considered which In this paper a forced, one-degree-of-freedom possesses both homoclinic and heteroclinic orbits (see section 2 for a description of the system, and section 3 for the definition of homoclinic and heteroclinic orbits). Classical approximate methods as well as the method of Melnikov have been used in the study. This study was motivated by the experimental work done by Moon [l] and the theoretical analysis done by Holmes [2] for a two-well potential dynamical system:
i+si-x+x3=
ycoswt.
The two-well potential system possesses homoclinic orbits. Holmes was able to calculate a necessary condition for the occurrence of strange attractor motions using the Melnikov method. In laboratory experiments, Moon showed that the homoclinic bifurcation criterion given by Holmes indeed gave a lower bound in the parameter space of the driving frequency and amplitude, (w, y). Using an approximate technique, Moon also derived a criterion in the (w, y) plane, which was in closer agreement with the experimental data. In this paper, the forced motion of an elastic beam with three stable and two unstable equilibrium positions, is a natural extension of the two-well potential system. A distinct 17 0022-460X/90/010017+ 18 %03.00/O
0 1990 Academic Press Limited
18
G. X.
LI AND
F. C.
MOON
feature of the system is that besides homoclinic orbits, it also posseses heteroclinic orbits. Thus one is able to study the system dynamics due to either homoclinic bifurcations or heteroclinic bifurcations or both. One of the purposes of this paper is to present a clear and fairly complete picture of the system chaos with relation to homoclinic and heteroclinic bifurcations. Each bifurcation will be shown to lead to a necessary chaos criterion. These results indicate that for the three-well problem, two criteria in the parameter space must be met for chaos in contrast to the two-well case where only one criterion is necessary. This study, and another by the first author in a doctoral dissertation involving a four-well potential problem, suggests that when there exists more than two attractors, more than one necessary condition may be required for chaotic vibrations. In section 2 we briefly describe the derivation of the equation of motion for an elastic beam with multiple equilibria. In the experimental system, this consisted of a cantilevered ferromagnetic beam with three magnets near the free end. The mathematical model of the experimental system must be accurate enough so that the predictions obtained from theoretical analysis make sense, and at the same time must also be simple enough to permit theoretical analysis. The stability properties of the system and the calculations of the homoclinic and heteroclinic bifurcations are performed in section 3. In section 4 numerical results are presented. The details of the numerical simulation of homoclinic and heteroclinic orbits and their bifurcations are discussed. To verify the correctness of the theoretically predicted bifurcations, the stable and unstable manifolds are also simulated so that one can observe when they intersect. Other features, such as Poincare maps, Lyapunov exponents and fractal basin boundaries are also simulated. In section 5 we present the experimental results for a chaos criterion in the plane of driving frequency and amplitude. Also in this section singular perturbation and harmonic balance methods are used to derive a functional relation between the forcing frequency and amplitude, with other parameters fixed, to predict the occurrence of chaos. We conclude our analyses in section 6.
2. EQUATION
OF MOTION
FOR THE EXPERIMENTAL
SYSTEM
The experimental system in this study consisted of an elastic ferromagnetic beam vertically clamped at the upper end to a stiff aluminum frame, and three permanent rare-earth magnets on the base plate of the frame (see Figure 1). The elastic beam has dimensions of 18.8 cm long, 9.5 mm wide and 0.23 mm thick. The whole system was shaken in the horizontal direction. The motion of the beam is considered to be planar. Under the magnetic forces generated by the magnets, the beam has three stable equilibrium positions with its lower-end tip near each magnet, and two unstable equilibrium positions between the stable ones. This problem is thus an analog of a particle in a three-well potential. The driving frequencies in our experiments are close to the lowest natural frequency of the system. The system dynamics is mainly dominated by the first mode of the beam. Therefore, we can formulate an ordinary differential equation to describe the dynamics of the system. When unforced, the beam vibrates in either one of the potential wells depending on initial conditions. Assuming that the damping is linear and proportional to the beam velocity, we can posit the governing differential equation in terms of the first mode of the beam, ij+&j(r)+crq(q2-b2)(q’-a*)=AR*cos(~t).
(2)
where q(f) represents the tip displacement of the beam, 8 the damping coefficient, a, b the stable and unstable equilibrium positions (b < a), l2 and A the driving frequency
THREE-WELL
POTENTIAL
OSCILLATOR
19
CHAOS
Buckled positions
Permar& Figure
1. Schematic
magnets
of the experimental
system.
and amplitude, and (Ya constant to be determined by the natural frequency of the beam at the buckled positions. Introducing the non-dimensional quantities x=9/a,
p = ab4/L’;,
s = ?7/n,,
one has the non-dimensional
x0= b/a,
w = n/n,,
f =A/a,
form of equation (2):
~+fs~+px(x~-x~)(x2-1)=f~*c0s#t.
(3)
For the purpose of the analytical study in the next section, this equation can be further resealed in the time variable and put in the form ~+~~+x(xz-x~)(x2-l)
= y cos ot,
(4)
where S is the damping ratio, o and y are equivalent to driving frequency and amplitude. Note that this is a Duffing’s type equation, non-linear up to fifth order in x. 3. THEORETICAL ANALYSIS PROPERTIES 3.1. STABILITY To investigate the stability properties of equation (4), we rewrite it in the first-order form in the phase plane
i = y,
p=-6y-x(x”-x:)(x2-l)+ycoswt.
(5)
To discuss the local and global stabilities, we first consider the unforced case: X?= y,
~=-sy-x(x2-x:)(x2-1).
(6)
This system is autonomous and has a phase portrait as shown in Figure 2 with and without damping. The local stabilities at the equilibrium positions can be easily determined by the eigenvalues of the linearized system. There exist three sinks, one at the origin with eigenvalues of A = (-8 f X/S’- 4x$/2( Re (A)
20
G. X.
LI AND
F. C.
MOON
Figure 2. Unforced phase flows of system (4); the two dotted closed orbits arc the homoclinic unperturbed system (6 = 0, y = 0), and the closed orbit in the middle is the heteroclinic orbit.
orbits for the
The information on global stability of the system can be obtained by looking for a Lyapunov function, which should be positive definite for a stable system. This is simple for the present system since the total energy of the system (kinetic and potential) is such a Lyapunov function: V=;y~+;x;x*-~(l+X;)X4+~X6.
(7)
Differentiating this function along the solution curve of system (6) gives d V/dt = -6y’, which is always negative for positive damping and y # 0. Therefore the global solutions are always bounded (see Figure 2). At the x-axis (y = 0), the tangents of the solutions are parallel to the y-axis since the x-component of the vector field is zero (a = 0) at such points. We also can expect that the forced system (5) is also globally bounded for small forcing y. To see this, let us once again differentiate the Lyapunov function (7) and this time along the solution curves of equation (5). The result is dV/dr=-6y’-6yycoswt~
-6(y’-ylyl).
It is clear that when forcing amplitude y is sufficiently small, or when y is sufficiently large, dV/dt can always be non-positive. In fact, for a moderate value of y the solution of system (5) will be either a periodic one or a non-periodic bounded one when t + +cc in the phase plane. 3.2. HOMOCLINIC AND HETEROCLINIC BIFURCATION CRITERIA FOR CHAOS The idea that homoclinic bifurcations in a dynamical system can be used to predict the occurrence of chaotic motions was originally introduced by Melnikov [3], and later was further developed by Holmes [2] and Salem [4]. The successful application of the method to the two-well potential dynamical system by Holmes [2] and Moon [l], and to the pendulum systems by Salem [4] and Moon et al. [S] show the usefulness of the method in non-linear dynamical systems. A homoclinic orbit in a plane is a closed curve intersecting at a saddle (hyperbolic) point. When an orbit connects two saddle points at its two ends then it is called a
THREE-WELL
POTENTIAL
OSCILLATOR
21
CHAOS
heteroclinic orbit. The three-well potential system (5), when unperturbed 8 = y = 0, possesses two homoclinic orbits and two heteroclinic orbits, as shown by dotted lines in Figure 2. We only need to study the dynamics of one homoclinic orbit and one heteroclinic orbit due to the system symmetry. When a T-periodic perturbation is added to the unperturbed system, then the homoclinic or heteroclinic orbits are referred to as the orbits in the Poincare map section (see section 4.2 for the details of numerically computing these orbits). Such kinds of orbits for the unperturbed system are structurally unstable: i.e., they can be destroyed by sufficiently small perturbations. The idea of bifurcations related to these orbits is that when the perturbations are added, it is possible that the broken orbits will reconnect again. Generally, such connections are transversal, and one transversal intersection will lead to infinitely many intersections in an extremely complicated way. If this happens in a system, chaos becomes possible [6,12]. Use of Poincare maps is one way to characterize the behavior of a chaotic motion in a two-dimensional phase plane. For a periodically forced two-dimensional dynamical system, like system (5), a natural choice of the Poincare section is one of the phase planes at a specified phase, {(x, y, 1)) t(mod T) = lo}. The Poincari map is the intersections of solutions of the forced system with the Poincare section. If the motion is periodic, then its Poincare map only show finitely many points. However, when the motion becomes chaotic, which is non-periodic but bounded, its Poincare map will contain infinitely many points. When a homoclinic or heteroclinic orbit breaks in a Poincare section, there are two stable manifolds, denoted by Ws, and two unstable manifolds, denoted by W”, across the corresponding saddle points (see Figure 4 of section 4.2). It is clear that when the perturbations are removed (6 = y = 0) the stable and unstable manifolds are identified, yielding either the homoclinic orbit or heteroclinic orbits. We point out that when motions of a dynamical system are bounded, the dynamics of the system due to homoclinic bifurcations is qualitatively the same as the dynamics due to heteroclinic bifurcations. Therefore we will treat the homoclinic orbits and heteroclinic orbits in our system qualitatively the same way. However, in this paper we want to emphasize the multiplicity of homoclinic or heteroclinic orbits: if two or more such orbits exist, what will happen when one or more than one or all bifurcations occur? We will use two dynamical concepts, Lyapunov exponents and fractal basin boundaries, to answer the above questions. We return to discuss the bifurcation problem. For a two-dimensional dynamical system, the differential equations can be put in the form a=fl(X,Y)+%l(%Y,
f, PL),
.+ =h(x,Y)+4X,Y,
t, PI,
(8a,
b)
where the unperturbed vector field (f,(x, y),f2(x, y)) is assumed to possess a homoclinic orbit in the phase plane, the perturbations (g,(x, y, t, p), g,(x, y, t, p)) are assumed to be T-periodic in time t, E is a sufficiently small number, and p is a free parameter. When the perturbations are added (E f O), the closed homoclinic orbit breaks. The central idea of the Melnikov method is to find a function which can be used to measure the separation between the stable and unstable manifolds, W’ and W”, along the unperturbed orbit. If this function vanishes for a certain parameter value of p then W” and W” will meet again away from the saddle point in the Poincare section. By a theorem attributed to Poincare (see reference [6]), if W‘ and W” cross each other once, they will intersect an infinite number of times. In such a case one can show that small regions in the Poincare section near the saddle point get stretched and folded in what is now called a horseshoe map. In a dissipative system, this stretching and folding leads to fractal properties of the motion in phase space and is a precursor for chaos. Suppose that the unperturbed homoclinic orbit can be found, (x0, yO)= (x,(t), y,Jr)), then it can be shown that the
G. X.
22
LI AND
F. C.
MOON
function +0? [fi(x,, A(&,) = I -K
y,)g,(xo, y,, t+ to, pL)-h(xo,
yo)g,(xo,
yo, t-t to, cL)I dt
(9)
does provide such a measure for sufficiently small E. (For more details, see reference [6].) Thus the vanishing of A in equation (9) becomes a necessary condition for chaotic vibrations. In order to apply the Melnikov formula to the dynamical system (5), we assume that 6 and y are small; specifically, we let 6 = E& y = ~7, where E c< 1 and d and i; are of order one. Thus system (5) can be put in the form i = y,
4;= -x(x’-
xi)( X2-l)+&(yCOS(wt)-6;).
(IO)
Comparing terms in expressions (10) and (S), one has f, =y, fi= -x(x’-x$(x’g, = 0, and g, = 7 cos wt - 8yo. Substituting into equation (9), gives +a +a3 A(t,) = -& ye(t) cos o(t+to) dt, y;(r) dt+ 7 I -m I -cc
l),
(11)
where to is the cross-section time of the Poincari map: to can also be interpreted as the initial time of the forcing term. The problem now becomes how.to compute yo( t), from system (5), for the homoclinic and heteroclinic orbits when 6 = y = 0. Unlike the two-well potential problem (l), where yO(t) can be found analytically, finding an analytical solution for yO(t) is impossible for the present three-well potential system. (Details of how to compute y,(t) numerically are given in the next section.) Note that y,(t) is a function of time from -co to +co. One must properly choose a starting point. With reference to Figure 2 again, a natural choice of the starting points would be p, for the homoclinic orbit and pZ for the heteroclinic orbit for the purpose of keeping symmetry. In this way, ye(t) would be an odd function of time for the homoclinic orbit and an even function for the heteroclinic orbit. Note that this will result in two necessary conditions or criteria for chaotic vibrations. We first compute the Melnikov function for the homoclinic orbit. Since yo( t) is odd in this case, equation (11) can be simplified to X uA( to) = -26 yi( t) dt - 27 sin wto yo( t) sin wt dt I0 I0 =
-2&
- 27 sin do Ihorn(
(12)
where a = 5: yi( t) dt is simply a constant once yo( t) is given, and Z,,om= 50”yo( t) sin wt dt is a function of the frequency w. Obviously, the necessary and sufficient condition for A ( to) to have simple zeros is 71,,, (w ) > 6a
or
Y’
‘Yhom
=
(13)
SalIhom(W).
For the heteroclinic
orbit, y,,(t) is even. Thus equation (11) can be reduced to ur A( to) = -26 y;(t) dt+27 cos wto y,,(t) cos wt, dt I0 I0 =
-2&Y +2-7 cos Wt&.,(W),
where Iher= jy yo( t) cos wt dt. Therefore, only if 7’
Yhet
the heteroclinic =
~~lbwr(W).
(14) bifurcations
will occur if and (15)
THREE-WELL
POTENTIAL
OSCILLATOR
4. COMPUTER
23
CHAOS
STUDIES
This section is purely for the purposes of numerical verification of the results obtained in the previous sections. Thus, we will choose the unstable equilibrium position x0 = 0.5 rather than 0.435 for the experimental system. The damping ratio S is chosen to be 0.1 (S = 0.03 in experiments). 4.1.
SIMULATION
OF
THE
BIFURCATION
CRITERIA
FOR
CHAOS
Our algorithms for computing the bifurcation conditions (13) and (15) basically include two steps. We first numerically solve the differential equation (6) for yO(r) using the fourth order Runge-Kutta technique, and then we numerically integrate the integrals (Y,Iho,,, and ZAP,. Even though the time intervals in those integrals are infinity, only finite time intervals are needed in our computations since yO(t)‘s exponentially converge to the corresponding saddle points, at which y,Jco) = 0. For our three-well potential system, the total length of time in the computations was 10 non-dimensional time units for both the homoclinic and heteroclinic orbits. The initial conditions at p1 and pz will give half of the homoclinic orbit and half of the heteroclinic orbit for forward time integrations. Once yO(I) is found, the integrals (Y, Zhomand Z,,e, are simply a matter of finite summations over the time interval O-10. The critical values of bifurcations are then obtained with the use of the formulas (13) and (15) for the given parameter values. Because the integration processes were involved, the time step 0.01 was used in the Runge-Kutta algorithm for computing yo(r)‘s. In other cases, the time step is taken to be 0.1. Note that the critical values for bifurcations, y,,Yhom and yher, are functions of frequency w. Saving the solutions of yO(t)‘s, one only needs to compute the integrals Zho,,,and I,,,,, and then YI,~,~and yher, for each value of w. For x0 = O-5 and S = 0.1, 50 homoclinic values and 50 heteroclinic values of y were computed over a frequency range of O-6-1.2. The results are shown in Figure 3, where the solid curve is for homoclinic bifurcations, while the dashed curve is for heteroclinic bifurcations. In the next few subsections, we will investigate what happens to the solutions as one crosses these bifurcation curves. 4.2. COMPUTATIONS OF THE STABLE AND UNSTABLE MANIFOLDS The computations of the homoclinic or heteroclinic orbits for the unperturbed system are straightforward. They are obtained by simply numerically integrating (6) with the proper starting points on the orbits. However, when the T-periodic perturbations are
o.ootl 0.6
0.7
0.9
0.9
I.0
1.1
I.2
Frequency
Figure
3. Homoclinic
(---)
and heteroclinic
(- - -) bifurcation
curves in the (w, y) plane;
x0 = 0.5, 6 = 0.1.
24
G.
X.
LI AND
F. C.
MOON
added, the computations of the manifolds have to be performed in the Poincare section. In this case the flow takes place in a three-dimensional phase space (x, y, wt) where the phase of the periodic force is taken modulo 27r/w. When the perturbations are sufficiently small, it can be shown that the Poincare map also has saddle points and stable and unstable manifolds which are close to the corresponding ones for the autonomous system (Y =O). To compute the manifolds in the Poincare section, one first has to find the corresponding saddle points. To be specific, suppose x0 is a saddle point for the following autonomous dynamical system i = f(x),
(16)
where x = (x, y) are the state variables, and f(x) is a vector field where f(xo) = 0. When a T-periodic forcing is added, i = f(x) + &g(t),
(17)
we expect that the saddle point will be slightly shifted for E <<1 in the Poincare section. To find the shifted saddle point x,, one solves the solution to (17) and then determines the PoincarC map by letting t = nT, n = 0, 1,2, . . . . In general the global Poincare map in the plane cannot be determined analytically. However, to find x,, one can linearize equation (17) in the neighbourhood of x0 by letting x = x0+ U: ti = Df(xo)u + &g(t) = Au + &g(t),
(18)
where A = Df(x,) is the Jacobian matrix of f(x) at x0 and is a constant matrix. By doing so, the solution to equation (18) can be obtained as I
to) + u(t) = e A(‘-‘,+L4(
E
I
eActp7)g( 7) dr.
(19)
‘0
One obtains the Poincare map in the neighbourhood u(T+to)=eAru(lo)+E
of x0 by letting t = T+ to
T+‘C, e A’T+‘dg(~) I ‘,I
d7,
and, thus, the perturbed saddle point x, =x,+e(Z-e
AT))’ I,““J eA(T+r,,-~)g(r) d7.
(20)
Note that x, is indeed E-close to the unperturbed one x0: 11x,-x0(/ = O(E). Once the perturbed saddle point x, is determined, one linearizes the vector field f(x) at x, and then determines the eigendirections of IIf( which give the tangents to the corresponding stable (with negative eigenvalues) and unstable (with positive eigenvalues) manifolds W” and W” crossing x,. To compute W”, we choose a small segment along the unstable eigendirection centred at x,. Many points along the segment are chosen to be initial conditions for equation (17), which can numerically integrated. The Poincari map is constructed from the solution of equation (17) by storing points which have time at nT (n=0,1,2,. . .). The resulting points lie on the unstable manifold of the Poincare map of the flow (17). The stable manifold W” can be computed in a similar way, with negative time steps. We now show some stable and unstable manifolds for system (5) in the Poincare map section (see Figures 4 and 5). We want to verify whether or not the stable manifold W” meets the unstable manifold in the Poincari map plane when the driving amplitude reaches the bifurcation value for which A = 0 in equation (12). When x0 =O*S, the
THREE-WELL
Figure 4. Homoclinic
bifurcations
POTENTIAL
OSCILLATOR
of system (4) in a Poincark
section;
CHAOS
x0 = 0.5,s
= 0. I, w = 1.O and y = 0.0334.
t Figure 5. Heteroclinic
bifurcations
of system (4) in a Poincark
section;
x0 = 0.5, 6 = 0.1, w = 10. and 7 = 0.12.
equilibrium point (x,, yO) = (0.5, O-O) is a saddle point for the unforced system (y = 0). For the fixed parameter values 8 = 0.1, w = 1.0, the saddle point for the Poincare map is moved to (0*476,0*00176) at the homoclinic bifurcation value y=O*O334, and to (0.4154,0-00615) at the heteroclinic bifurcation value y =O-117. As shown in Figure 4 for the homoclinic bifurcations, W’ just meets W” tangentially at a few points at the predicted critical value of y = OeO334.For the heteroclinic bifurcations (see Figure 5), it appears that when y = O-12, w” and W” meet, compared with the predicted value of O-117. Note that while the homoclinic bifurcation occurs at y = O-0334 in Figure 4, the stable and unstable manifolds, corresponding to the heteroclinic orbits, are still far away from each other.
G. X.
26 4.3.
LYAPUNOV
LI
AND
F. C.
MOON
EXPONENTS
Having obtained the critical values of the homoclinic and heteroclinic bifurcations, we are in a position to search for the system chaos. We want to investigate, if a homociinic or heteroclinic bifurcation occurs in the system, what will happen with the system response. We answer this question by employing the Lyapunov exponent technique in the driving frequency and amplitude plane (see, e.g. reference [12]). If a motion is chaotic, then all neighboring motions will diverge from it exponentially. The Lyapunov exponent is simply a constant measuring such divergence or convergence on average. For a two-dimensional dynamical system, like system (5), there are two Lypunov exponents associated with each solution. If at least one exponent is positive, then the motion is said to be chaotic. For our purpose of finding chaos, only the largest exponent is needed and simulated for each motion. The advantage of the Lyapunov exponent technique in the study of chaos is that one can systematically investigate the chaotic behavior of the system in a region in a parameter space by simulating the motion and its Lyapunov exponent simultaneously. For the dynamical system (5), we calculated 100 x 100 Lyapunov exponents in the (w, y) plane for the initial condition (0.5,O.l). Only those parameter values of (w, y) with positive exponents were stored and plotted. In Figure 6, each + represents a positive exponent with value greater than, or equal to 0.01 at the given driving frequency and amplitude. The simulation result immediately tells us that many, but not all, parameter values above these bifurcation curves do cause chaotic responses in the system. The higher the driving amplitude is, the more likely chaos is found. We also carried out similar simulations for the initial condition (0.49,0.0), and the result was found to be qualitatively the same as the one in Figure 6.
++
J+
+:+
+
+
+
&,A
i’+*jd+++++ / + + ++ + + +
+ 0.08 -
/
/
+-
/
/
/
/
+,
++/ j4+ t + + + :+$+ / ++
+ t
/’
0.06 -
*++1
//+ ++,#
+
_
+ +
+
/’
+ _ +
_
, /’ 0.04 >./’ _ + 0.02 -
Figure 6. Positive Lyapunov bifurcation curves: 8 = 0.1.
t
.H’ +
exponents
in the (w. y) plane;
chaotic
motions
become
possible
above
the
The Lyapunov exponents in Figure 6 were simulated in a Micro VAX computer. The integration step size was 0.1, and each exponent was averaged in a time period of 650 non-dimensional time units. By checking some of those positive exponents marked by + in Figure 6, we found that they were indeed positive within the time period. However, if the averaging time period is larger than 650, no guarantee can be given to their final values. In fact, experience tells us that an apparent chaotic motion sometimes escapes to a regular pattern after a sufficiently long term. In the next subsection, we present two
THREE-WELL
POTENTIAL
OSCILLATOR
27
CHAOS
phase portraits for some special parameter values in Figure 6 to see the phenomenon transient chaos and persistent chaos. 4.4.
PHASE
of
TRAJECTORIES
To see the effect of the homoclinic and heteroclinic bifurcations, several parameter values of (w, y) were chosen across the bifurcation curves in Figure 6, and the phase portraits were simulated for damping 6 =O*l and the initial condition of (x0, y(,) = (0.5,O.l). For (w, y) = (O-918,0*0435), which is just above the homoclinic bifurcation value of (0.918,0*0312), the phase trajectory is shown in Figure 7(a). The time period for the simulation was 600 non-dimensional time units. We found that the motion appeared random-like in the first 200 time units and then gradually converged to a limit cycle. We also monitored the Lyapunov exponent simulation together with the phase trajectory, and found that the exponent converged approximately to 0.12 for the first 200 time units and then slowly decayed toward zero. In the literature, such a motion is referred to as transient
chaos.
To find a permanent chaos in the time period of 600 time units, we chose another point, (O-714,0*132) from the middle of the cloud in the (w, y) plane (see Figure 6), for simulating another phase protrait. In Figure 7(b) is shown the phase trajectory for this case, and one can clearly see that the motion was indeed chaotic within this time period. Therefore, for our three-well potential system, numerical evidence in Figure 7 suggests that both the homoclinic and heteroclinic bifurcations are necessary conditions for the
r. 5 9
0.4
0.6
0.5
0.7
0.6
-0.6 -1.5
-1.0
-0.5
0.0
0.9
I.0
0.5
1.1
I.0
1.2
I.5
Disdacement
Figure 7. Phase trajectories; (a) transient chaos, 6 = 0.1, w = 0.918, y = 0.0435; (b) persistent chaos, S = 0.1, w =0.714, y =0.132. The simulations were carried out in a period of 600 time units with initial condition at (0.5,O.l) for both cases.
28
G. X. LI AND
F. C. MOON
occurrence of long-term chaotic motions of the system. Satisfaction of conditions for at least one of these bifurcation appears necessary for transient chaos to occur. 4.5.
FRACTAL
BASIN
BOUNDARIES
Basin boundary simulations are another means for checking homoclinic or heteroclinic bifurcations [9, lo]. When a homoclinic bifurcation occurs, the stable manifold W” will (a)
Figure
8. Basin boundaries:
(a) smooth,
y =0.02;
(b) fractal,
y =0.034;
(c) more fractal,
y=O.O5;
S =O.l.
THREE-WELL
POTENTIAL
OSCILLATOR
CHAOS
29
intersect with the unstable manifold W” infinitely many times in the Poincare section. Suppose that the dynamical system has two or more regular attractors, such as two limit cycles for example. A choice of two initial conditions on either side of a stable manifold will lead to two different final states of motions. In fact, the basin boundary in the initial condition plane is the closure of the stable manifolds of the system, just like the Poincare map is the closure of the unstable manifolds. Homoclinic or heteroclinic bifurcations lead to complicated intersection of the stable and unstable manifolds. Therefore, at these bifurcations the basin boundaries will become extremely complicated too, and are called fructal basin boundaries. Such fractal boundaries indicate that whether the system is attracted to one or or the other periodic attractor may be very sensitive to initial conditions. Thus small uncertainties in the initial conditions can lead to unpredictability of the system output even when the motion is not chaotic. This sensitivity to initial conditions, however, may become a precursor to chaotic dynamics. Numerical simulation was carried out in order to determine the basin boundaries for the three-well potential system (5) for different values of driving amplitudes y at S = 0.1, o = 1*O(see Figure 8). The system has three attractors in the neighborhoods of x = 0, * 1. For small y, the steady state motion will be a periodic limit cycle in one of the three attractors. Even if y is increased beyond the critical values of bifurcations, it is still possible that the final steady motion could be periodic rather than chaotic. In fact, near the bifurcation curves (see Figure 6), no permanent chaos was found in our simulations. The colors blue, yellow, and red were used in our simulations to denote initial conditions (starting points), the final steady states of which converged to the left (x = l), middle (x = 0), and right (x = 1) attractors. In Figure 8 are the black and white versions of the original color simulations. In Figure 8(a), y =0.02, which is smaller than the critical value 0.0334 for the homoclinic bifurcations. The boundaries are clearly smooth. Figure 8(b) corresponds to the predicted critical value 0.0334 (= yhom) at w = 1.0. Here many sections of the basin boundaries become non-smooth or fractal. If initial conditions are chosen near the fractal basin boundaries with a small uncertainty in this case, then the final states to which the motion is attracted is not predictable. The uncertainty of final states in the presence of small variations in the initial conditions is a feature of chaotic systems. With further increases in y, the fractal basin boundaries spread out and become more complex, as shown in Figure 8(c) for y = 0.05. By now it is clear that a homoclinic or heteroclinic bifurcation is a necessary condition for chaotic responses of our three-well potential system and a sufficient condition for the appearance of fractal basin boundaries. Numerical simulations tell us that the probability of finding chaos is higher for reasonably larger driving amplitudes.
5. EXPERIMENTAL CHAOS OF THE NONLINEAR SYSTEM In this section we investigate the chaotic responses for the experimental system. Our purposes are twofold: (1) to see if chaos exists in the system; and (2) if it exists, to see if the mathematical model (5) can be used to predict chaos. The external forcing of the experimental system was sinusoidal, supplied by an electric shaker. To measure the motion of the beam, two strain gages were mounted on the upper end of the beam near the clamped end, which produced electrical signals proportional to the bending of the beam. The velocity signals were obtained with the use of an electronic differentiator. The experimental set-up is shown in Figure 9. Without the magnets, the lowest three natural frequencies of the beam were measured to be 4.6, 26.6 and 73.6 Hz. When the magnets were added, the beam had three stable equilibrium positions. The lowest natural frequencies of the beam at the three potential
30
C.
Figure 9. Schematic
X.
LI
AND
F. C. MOON
of the data acquisition
system
of experiments.
wells were found to be 7.8, 8.1 and 8.3 Hz in the left, middle, and right potential wells respectively; they were not symmetric as expected. The asymmetry might be due to prestress in the beam and/or different magnetic strength in the magnets. The experiments were performed at the driving frequencies from 5 to 10 Hz, which were close to the three lowest natural frequencies of the beam. When the driving amplitude was small for a fixed frequency, the response of the beam was periodic of period one. However, when the driving amplitude was increased, we often observed subharmonic motions with period two or more. When the amplitude was increased beyond some critical value, the response would become chaotic. In our experimental studies, Poincare maps were used as a means to distinguish chaotic signals from periodic motions. If a motion was periodic, the Poincare map in the phase plane contained a finite number of points. When the motion was chaotic, there appeared a fractal-looking Poincare map in the (x, li) plane. This highly ordered, fractal-like structure is a basic feature distinguishing chaos from random noise, and is often called a strange attractor. The procedure to experimentally construct a Poincare map can be found in a book by the second author [12]. For the three-well potential system, a typical Poincare map for a chaotic motion is shown in Figure 10. For the purpose of comparison, a numerical Poincare map was also simulated by using equation (3). The parameter values used were /3 =0.577, x0=0.435, w =0.833, which correspond to the experiment, and 6 = 0.15 and f = 0.36, which were higher than the measured ones of S = 0.03 and f= 0.08. We should point out that it was difficult to generate numerically a Poincare map for small damping 6. When 6 is small, it is difficult to find a persistent chaotic response (long-time simulation) which is needed to generate a Poincare map with sufficiently many points. Also, when 6 is small, the Poincare map is very diffuse (close to a forced Hamiltonian system) and hardly any fractal pattern can be observed. In Figure 11, 4000 points (corresponding to 4000 forcing periods) were generated. To find regions in the driving frequency-amplitude plane in which motions of the dynamical system are chaotic, the driving amplitude was slowly increased for a fixed driving frequency near the natural frequencies. When the Poincari map of a response became fractal-looking, like the one in Figure 10, then the motion was declared to be
THREE-WELL
Figure oscillator
10. An experimental is 0.03.
PoincarC
POTENTIAL
map of chaotic
motion;
31
CHAOS
OSCILLATOR
the damping
ratio measured
from the physical
0.6-
2, 0.4.E s 9
0.2-
o.o-
-0.2-
-0.4~""""""""""""""~7 -1.5 -1.0 -0.5
0.0
0.5
I.0
1.5
Dlsplocement Figure 11. A PoincarC w = 0.833, j’= 0.36.
map
from
numerical
simulation;
in equation
(3), 6 =0.15,
/3 =0.577,
x,=0.435,
chaotic at that particular parameter value. Repeating the above procedure for different driving frequencies in the range of 5-10 Hz, we obtained an experimental threshold for chaos in the driving frequency and amplitude plane (see Figure 12). We point out that all motions in the experiments had initial conditions in the left potential well. The duration of observation of a chaotic motion was about four minutes. We also should point out that no guarantee can be given whether or not those chaotic motions are still chaotic at a later time. Also, if a response of the system is chaotic for a particular initial condition, it may not be chaotic for other initial conditions. In Figure 12, we also plot the corresponding homoclinic (dashed curve) and heteroclinic (solid curve) bifurcations by converting the non-dimensional quantities into their physical counterparts. The homoclinic bifurcation curve is considerably lower than the experimentally measured one. The heteroclinic one seems much better, especially at higher frequency. It is clear that while the analytical criteria (13) and (15) are satisfied, these
32
G.
X.
LI
AND
F.
C.
MOON
121,,,,[,,,,,,(,,
,(,,
,,,,
i
\
----o~-~~~.t..~.-c~.5
---__
__
_ -II--%__ 6
7
6
Drlvmg
-----_ I 9
IO
frequency
Figure 12. Criteria for chaos for the three-well potential dynamical formula (22); -, heteroclinic bifurcation (15); - - -, homoclinic
system. bifurcation
-, Experimental;
-
-, heuristic
(13).
lower bounds lie too far below the experimental results. In the following, we derive another criterion from a different point of view. Following an idea introduced in the study of the two-wall potential system by the second author [2], we also calculated a similar heuristic criterion to predict chaos for the three-well potential oscillator. The basic idea is to look for when periodic orbits become large enough to jump out of one of the potential wells (see also reference [12]). As mentioned earlier, the response of the system is periodic for small external forcing. A relation between the driving frequency and amplitude, also involving the response amplitude, could be found through the periodic-solution-solving process for small forcing. Based on experimental observations, we found that when the phase plane orbit become large enough, about the size of the unforced separatrices (the homoclinic and heteroclinic orbits in Figure 2), the motion often had a strange behavior leading to chaos. This property of the solutions was used to obtain a criterion for predicting chaos similar to the homoclinic/heteroclinic criteria (13) and (15). Since all chaotic motions in our experiments had initial conditions in the left potential well, we rewrite equation (3) about the left potential minimum, x = -1. This can be accomplished by replacing x with x - 1, x+&i+
f
;=,
aix’=fw2COS(Wt+#J),
(21)
where (Yr=2P(l
-xi),
a2 = -/?(7 -3x&
a3
=
P (9
-
d,,
a4
=
-sp,
(Yg=p.
Here a phase difference 4 has been introduced for the convenience of mathematical manipulations. In the master’s thesis of the first author [ 111, the harmonic method and the Lindstedt singular perturbation method have been used for the purpose to solve for an approximate periodic solution of the system under the assumption that f is small and w is close to the non-dimensional natural frequency. Doing so, we obtained f=(a/~~){~~~~+[&~~+(~(~~-5c&6cz,)u*]~}”~
(22a)
from the harmonic balance method, and (22b)
THREE-WELL
POTENTIAL
OSCILLATOR
CHAOS
33
from the singular perturbation method. In equations (22), w0 = J2p( 1 -x$) is the normalized natural frequency at the left potential well, and a is the amplitude of the fundamental component of the system response. Notice that both of equations (22) were derived under the assumption that 1w - q,I cc 1, and in equation (22b) we assumed x = a,+ a cos wt + a, cos 2wl, where a, 2 ---~cY~u~/cY,and a, = $x2a2/cx, following the analysis of the relative orders of magnitude. Obviously, under this assumption, the two equations are almost identical. When w is tuned away from wO, the error will be expected. The reason for deriving two relations at the same time is that when w is not close to w,,, probably each of the two relations is good for certain frequency ranges. This is the case in our analysis; the driving frequency can be quite far away from the natural frequencies. This criterion is based on the heuristic idea that chaos might occur when the velocity of the forced periodic motion attains a critical value. If a is a measure of the amplitude of motion, then the critical velocity is chosen to make the orbit in the phase plane close to the size of the separatrix of the unforced, undamped motion. This criterion proved successful in the case of the two-well potential problem (see reference [I]). To represent the closeness of the response to the separatrices, we use the velocity of the beam as a parameter, which can be approximated by aw(=u). If we once again use a “cheat” factor (Y(O< (Y< l), the above idea can be qualified by U=(YU”
or
wa = au0
(23)
near the chaotic regions, where z10is the maximum velocity along the separatrices for the unperturbed system. Substituting it into equations (22), we obtain relations between w and f, with LYas an adjustable constant to match with the experimental data. We found that by properly choosing the factor (~(=0.87), the relation (22a) compares well with the experimental data for w < wO, while the relation (22b) compares well with the data for w > wg. In Figure 12, these two curves are matched at the natural frequency w0 and give a much better prediction to the region of chaotic motions.
6. CONCLUSIONS that homoclinic or The results in this paper, once again, lead to the conclusion heteroclinic bifurcations in a dynamical system can be accurately predicted with the use of Melnikov’s method (see Figures 4 and 5). Chaotic motions of the three-well potential system were found both experimentally and numerically. A horseshoe map is found for the system when homoclinic or heteroclinic bifurcations occur, which often causes complex dynamics. It is important to be aware that although the existence of horseshoes gives the possibility of chaos, one cannot conclude the existence of chaos without additional information. The Poincare map technique is the best means to check the existence of a strange attractor motion in a laboratory. If the motion is really chaotic, then a highly ordered, fractal-looking, pattern should be observed. On the other hand, the Lyapunov exponent technique is very useful in checking chaos in numerical simulations. When doing so, care must be taken regarding the convergence of the exponent. The three-well potential system has three limit cycle attractors for sufficiently small forcing. The existence of homoclinic or heteroclinic bifurcations implies the existence of fractal basin boundaries of attracting domains and the loss of absolute predictability in the face of small uncertainties in initial conditions. In fact, the Melnikov technique is accurate in predicting the occurrence of such fractal basin boundaries. The experiments of the three-well potential system also showed chaotic motions, even when the forcing inputs were deterministic. In the frequency-amplitude plane, the experimental criterion for chaos is above both the homoclinic and heteroclinic bifurcation
G. X. LI
34
AND
F. C. MOON
curves, indicating that these bifurcations are necessary conditions for the appearance of chaos in the system. The heuristic criterion (22), however, gives a much closer prediction
for chaos.
REFERENCES
of Applied Mechanics 47, 638-644. Experiments on chaotic motions of a forced nonlinear oscillator: strange attractor. 2. P. J. HOLMES 1979 Philosophical Transactions of the Royal Society A (London) 292, 419-448. 1. F. C. MOON 1980 Journa/
A nonlinear oscillator with a strange attractor. 3. V. K. MELNIKOV 1963 Transacrions of the Moscow Mathematical
4. 5. 6. 7. 8. 9. 10.
Sociefy 12, l-57. On the stability of the center for time periodic motion. F. C. SALAM 1986 SIAM Journal of Applied Mathematics 47 (2), 232-243. The Melnikov technique for highly dissipative systems. F. C. MOON, J. CUSUMANO and P. J. HOLMES 1987 Physica 24D, 383-390. Evidence for homoclinic orbits as a precursor to chaos in a magnetic pendulum. J. GUCKENHEIMER and P. HOLMES 1983 Nonlinear Oscillafions, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag. A. WOLF, J. Swrrr, H. SWINNEY and J. VASTANO 1985 Physica 16D, 285-317. Determining Lyapunov exponents from a time series. A. J. LICHTENBERG and M. A. LIHERMAN 1983 Regular and Stochastic Motion. New York: Springer-Verlag. F. C. MOON and G. X. LI 1985 Physica/ Review Letters 55, 1439-1442. Fractal basin boundaries and homoclinic orbits for periodic motion in a two-well potential. C. GREBOGI, S. MACDONALI), E. OTT and J. YORKE 1983 Phvsics Letters 99A, 415-418.
Final state sensitivity: an observation to predictability. 11. G. X. LI 1984 M. S. Thesis, Cornell University, Ithaca, N. Y. Chaotic System with Five Equilibrium States. 12. F. C. MOON 1987 Chaotic Vibrations. New York: John Wiley.
Vibrations
of a Nonlinear
of Sound and
Vibration (1990) 136(l), 17-34
CRITERIA POTENTIAL
FOR
CHAOS
OSCILLATOR
AND
OF A THREE-WELL WITH
HETEROCLINIC G.
HOMOCLINIC ORBITS
X. Lr
Department of Theoretical and Applied Mechanics AND
F. C.
MOON
Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, New York 14853, U.S.A. (Received 13 July 1988, and in revised form 10 April 1989)
Criteria for chaotic vibrations of a non-linear elastic beam having three stabie and two unstable equilibrium positions are compared. The responses of the system exhibit chaotic behavior in certain regions in the driving frequency and amplitude plane. With the use of the Melnikov technique, necessary conditions for chaos are derived based on homoclinic and heteroclinic bifurcations of the system. Numerical simulations of the stable and unstable manifolds, Poincare maps, Lyapunov exponents and fractal basin boundaries are performed. The results indicate that both homoclinic and heteroclinic bifurcations are necessary conditions for chaos; thus, at least two necessary conditions must be satisfied for any chaotic responses. However, the final steady state of motion for many initial conditions becomes unpredictable (in the sense of fractal basin boundaries) if at least one homoclinic or heteroclinic bifurcation occurs. Finally, an experimental criterion for chaotic motion, and an heuristic criterion based on classical methods in non-linear dynamics are also obtained and compared with the bifurcation criteria, as a function of forcing amplitude and frequency.
1. INTRODUCTION
non-linear oscillator is considered which In this paper a forced, one-degree-of-freedom possesses both homoclinic and heteroclinic orbits (see section 2 for a description of the system, and section 3 for the definition of homoclinic and heteroclinic orbits). Classical approximate methods as well as the method of Melnikov have been used in the study. This study was motivated by the experimental work done by Moon [l] and the theoretical analysis done by Holmes [2] for a two-well potential dynamical system:
i+si-x+x3=
ycoswt.
The two-well potential system possesses homoclinic orbits. Holmes was able to calculate a necessary condition for the occurrence of strange attractor motions using the Melnikov method. In laboratory experiments, Moon showed that the homoclinic bifurcation criterion given by Holmes indeed gave a lower bound in the parameter space of the driving frequency and amplitude, (w, y). Using an approximate technique, Moon also derived a criterion in the (w, y) plane, which was in closer agreement with the experimental data. In this paper, the forced motion of an elastic beam with three stable and two unstable equilibrium positions, is a natural extension of the two-well potential system. A distinct 17 0022-460X/90/010017+ 18 %03.00/O
0 1990 Academic Press Limited
18
G. X.
LI AND
F. C.
MOON
feature of the system is that besides homoclinic orbits, it also posseses heteroclinic orbits. Thus one is able to study the system dynamics due to either homoclinic bifurcations or heteroclinic bifurcations or both. One of the purposes of this paper is to present a clear and fairly complete picture of the system chaos with relation to homoclinic and heteroclinic bifurcations. Each bifurcation will be shown to lead to a necessary chaos criterion. These results indicate that for the three-well problem, two criteria in the parameter space must be met for chaos in contrast to the two-well case where only one criterion is necessary. This study, and another by the first author in a doctoral dissertation involving a four-well potential problem, suggests that when there exists more than two attractors, more than one necessary condition may be required for chaotic vibrations. In section 2 we briefly describe the derivation of the equation of motion for an elastic beam with multiple equilibria. In the experimental system, this consisted of a cantilevered ferromagnetic beam with three magnets near the free end. The mathematical model of the experimental system must be accurate enough so that the predictions obtained from theoretical analysis make sense, and at the same time must also be simple enough to permit theoretical analysis. The stability properties of the system and the calculations of the homoclinic and heteroclinic bifurcations are performed in section 3. In section 4 numerical results are presented. The details of the numerical simulation of homoclinic and heteroclinic orbits and their bifurcations are discussed. To verify the correctness of the theoretically predicted bifurcations, the stable and unstable manifolds are also simulated so that one can observe when they intersect. Other features, such as Poincare maps, Lyapunov exponents and fractal basin boundaries are also simulated. In section 5 we present the experimental results for a chaos criterion in the plane of driving frequency and amplitude. Also in this section singular perturbation and harmonic balance methods are used to derive a functional relation between the forcing frequency and amplitude, with other parameters fixed, to predict the occurrence of chaos. We conclude our analyses in section 6.
2. EQUATION
OF MOTION
FOR THE EXPERIMENTAL
SYSTEM
The experimental system in this study consisted of an elastic ferromagnetic beam vertically clamped at the upper end to a stiff aluminum frame, and three permanent rare-earth magnets on the base plate of the frame (see Figure 1). The elastic beam has dimensions of 18.8 cm long, 9.5 mm wide and 0.23 mm thick. The whole system was shaken in the horizontal direction. The motion of the beam is considered to be planar. Under the magnetic forces generated by the magnets, the beam has three stable equilibrium positions with its lower-end tip near each magnet, and two unstable equilibrium positions between the stable ones. This problem is thus an analog of a particle in a three-well potential. The driving frequencies in our experiments are close to the lowest natural frequency of the system. The system dynamics is mainly dominated by the first mode of the beam. Therefore, we can formulate an ordinary differential equation to describe the dynamics of the system. When unforced, the beam vibrates in either one of the potential wells depending on initial conditions. Assuming that the damping is linear and proportional to the beam velocity, we can posit the governing differential equation in terms of the first mode of the beam, ij+&j(r)+crq(q2-b2)(q’-a*)=AR*cos(~t).
(2)
where q(f) represents the tip displacement of the beam, 8 the damping coefficient, a, b the stable and unstable equilibrium positions (b < a), l2 and A the driving frequency
THREE-WELL
POTENTIAL
OSCILLATOR
19
CHAOS
Buckled positions
Permar& Figure
1. Schematic
magnets
of the experimental
system.
and amplitude, and (Ya constant to be determined by the natural frequency of the beam at the buckled positions. Introducing the non-dimensional quantities x=9/a,
p = ab4/L’;,
s = ?7/n,,
one has the non-dimensional
x0= b/a,
w = n/n,,
f =A/a,
form of equation (2):
~+fs~+px(x~-x~)(x2-1)=f~*c0s#t.
(3)
For the purpose of the analytical study in the next section, this equation can be further resealed in the time variable and put in the form ~+~~+x(xz-x~)(x2-l)
= y cos ot,
(4)
where S is the damping ratio, o and y are equivalent to driving frequency and amplitude. Note that this is a Duffing’s type equation, non-linear up to fifth order in x. 3. THEORETICAL ANALYSIS PROPERTIES 3.1. STABILITY To investigate the stability properties of equation (4), we rewrite it in the first-order form in the phase plane
i = y,
p=-6y-x(x”-x:)(x2-l)+ycoswt.
(5)
To discuss the local and global stabilities, we first consider the unforced case: X?= y,
~=-sy-x(x2-x:)(x2-1).
(6)
This system is autonomous and has a phase portrait as shown in Figure 2 with and without damping. The local stabilities at the equilibrium positions can be easily determined by the eigenvalues of the linearized system. There exist three sinks, one at the origin with eigenvalues of A = (-8 f X/S’- 4x$/2( Re (A)
20
G. X.
LI AND
F. C.
MOON
Figure 2. Unforced phase flows of system (4); the two dotted closed orbits arc the homoclinic unperturbed system (6 = 0, y = 0), and the closed orbit in the middle is the heteroclinic orbit.
orbits for the
The information on global stability of the system can be obtained by looking for a Lyapunov function, which should be positive definite for a stable system. This is simple for the present system since the total energy of the system (kinetic and potential) is such a Lyapunov function: V=;y~+;x;x*-~(l+X;)X4+~X6.
(7)
Differentiating this function along the solution curve of system (6) gives d V/dt = -6y’, which is always negative for positive damping and y # 0. Therefore the global solutions are always bounded (see Figure 2). At the x-axis (y = 0), the tangents of the solutions are parallel to the y-axis since the x-component of the vector field is zero (a = 0) at such points. We also can expect that the forced system (5) is also globally bounded for small forcing y. To see this, let us once again differentiate the Lyapunov function (7) and this time along the solution curves of equation (5). The result is dV/dr=-6y’-6yycoswt~
-6(y’-ylyl).
It is clear that when forcing amplitude y is sufficiently small, or when y is sufficiently large, dV/dt can always be non-positive. In fact, for a moderate value of y the solution of system (5) will be either a periodic one or a non-periodic bounded one when t + +cc in the phase plane. 3.2. HOMOCLINIC AND HETEROCLINIC BIFURCATION CRITERIA FOR CHAOS The idea that homoclinic bifurcations in a dynamical system can be used to predict the occurrence of chaotic motions was originally introduced by Melnikov [3], and later was further developed by Holmes [2] and Salem [4]. The successful application of the method to the two-well potential dynamical system by Holmes [2] and Moon [l], and to the pendulum systems by Salem [4] and Moon et al. [S] show the usefulness of the method in non-linear dynamical systems. A homoclinic orbit in a plane is a closed curve intersecting at a saddle (hyperbolic) point. When an orbit connects two saddle points at its two ends then it is called a
THREE-WELL
POTENTIAL
OSCILLATOR
21
CHAOS
heteroclinic orbit. The three-well potential system (5), when unperturbed 8 = y = 0, possesses two homoclinic orbits and two heteroclinic orbits, as shown by dotted lines in Figure 2. We only need to study the dynamics of one homoclinic orbit and one heteroclinic orbit due to the system symmetry. When a T-periodic perturbation is added to the unperturbed system, then the homoclinic or heteroclinic orbits are referred to as the orbits in the Poincare map section (see section 4.2 for the details of numerically computing these orbits). Such kinds of orbits for the unperturbed system are structurally unstable: i.e., they can be destroyed by sufficiently small perturbations. The idea of bifurcations related to these orbits is that when the perturbations are added, it is possible that the broken orbits will reconnect again. Generally, such connections are transversal, and one transversal intersection will lead to infinitely many intersections in an extremely complicated way. If this happens in a system, chaos becomes possible [6,12]. Use of Poincare maps is one way to characterize the behavior of a chaotic motion in a two-dimensional phase plane. For a periodically forced two-dimensional dynamical system, like system (5), a natural choice of the Poincare section is one of the phase planes at a specified phase, {(x, y, 1)) t(mod T) = lo}. The Poincari map is the intersections of solutions of the forced system with the Poincare section. If the motion is periodic, then its Poincare map only show finitely many points. However, when the motion becomes chaotic, which is non-periodic but bounded, its Poincare map will contain infinitely many points. When a homoclinic or heteroclinic orbit breaks in a Poincare section, there are two stable manifolds, denoted by Ws, and two unstable manifolds, denoted by W”, across the corresponding saddle points (see Figure 4 of section 4.2). It is clear that when the perturbations are removed (6 = y = 0) the stable and unstable manifolds are identified, yielding either the homoclinic orbit or heteroclinic orbits. We point out that when motions of a dynamical system are bounded, the dynamics of the system due to homoclinic bifurcations is qualitatively the same as the dynamics due to heteroclinic bifurcations. Therefore we will treat the homoclinic orbits and heteroclinic orbits in our system qualitatively the same way. However, in this paper we want to emphasize the multiplicity of homoclinic or heteroclinic orbits: if two or more such orbits exist, what will happen when one or more than one or all bifurcations occur? We will use two dynamical concepts, Lyapunov exponents and fractal basin boundaries, to answer the above questions. We return to discuss the bifurcation problem. For a two-dimensional dynamical system, the differential equations can be put in the form a=fl(X,Y)+%l(%Y,
f, PL),
.+ =h(x,Y)+4X,Y,
t, PI,
(8a,
b)
where the unperturbed vector field (f,(x, y),f2(x, y)) is assumed to possess a homoclinic orbit in the phase plane, the perturbations (g,(x, y, t, p), g,(x, y, t, p)) are assumed to be T-periodic in time t, E is a sufficiently small number, and p is a free parameter. When the perturbations are added (E f O), the closed homoclinic orbit breaks. The central idea of the Melnikov method is to find a function which can be used to measure the separation between the stable and unstable manifolds, W’ and W”, along the unperturbed orbit. If this function vanishes for a certain parameter value of p then W” and W” will meet again away from the saddle point in the Poincare section. By a theorem attributed to Poincare (see reference [6]), if W‘ and W” cross each other once, they will intersect an infinite number of times. In such a case one can show that small regions in the Poincare section near the saddle point get stretched and folded in what is now called a horseshoe map. In a dissipative system, this stretching and folding leads to fractal properties of the motion in phase space and is a precursor for chaos. Suppose that the unperturbed homoclinic orbit can be found, (x0, yO)= (x,(t), y,Jr)), then it can be shown that the
G. X.
22
LI AND
F. C.
MOON
function +0? [fi(x,, A(&,) = I -K
y,)g,(xo, y,, t+ to, pL)-h(xo,
yo)g,(xo,
yo, t-t to, cL)I dt
(9)
does provide such a measure for sufficiently small E. (For more details, see reference [6].) Thus the vanishing of A in equation (9) becomes a necessary condition for chaotic vibrations. In order to apply the Melnikov formula to the dynamical system (5), we assume that 6 and y are small; specifically, we let 6 = E& y = ~7, where E c< 1 and d and i; are of order one. Thus system (5) can be put in the form i = y,
4;= -x(x’-
xi)( X2-l)+&(yCOS(wt)-6;).
(IO)
Comparing terms in expressions (10) and (S), one has f, =y, fi= -x(x’-x$(x’g, = 0, and g, = 7 cos wt - 8yo. Substituting into equation (9), gives +a +a3 A(t,) = -& ye(t) cos o(t+to) dt, y;(r) dt+ 7 I -m I -cc
l),
(11)
where to is the cross-section time of the Poincari map: to can also be interpreted as the initial time of the forcing term. The problem now becomes how.to compute yo( t), from system (5), for the homoclinic and heteroclinic orbits when 6 = y = 0. Unlike the two-well potential problem (l), where yO(t) can be found analytically, finding an analytical solution for yO(t) is impossible for the present three-well potential system. (Details of how to compute y,(t) numerically are given in the next section.) Note that y,(t) is a function of time from -co to +co. One must properly choose a starting point. With reference to Figure 2 again, a natural choice of the starting points would be p, for the homoclinic orbit and pZ for the heteroclinic orbit for the purpose of keeping symmetry. In this way, ye(t) would be an odd function of time for the homoclinic orbit and an even function for the heteroclinic orbit. Note that this will result in two necessary conditions or criteria for chaotic vibrations. We first compute the Melnikov function for the homoclinic orbit. Since yo( t) is odd in this case, equation (11) can be simplified to X uA( to) = -26 yi( t) dt - 27 sin wto yo( t) sin wt dt I0 I0 =
-2&
- 27 sin do Ihorn(
(12)
where a = 5: yi( t) dt is simply a constant once yo( t) is given, and Z,,om= 50”yo( t) sin wt dt is a function of the frequency w. Obviously, the necessary and sufficient condition for A ( to) to have simple zeros is 71,,, (w ) > 6a
or
Y’
‘Yhom
=
(13)
SalIhom(W).
For the heteroclinic
orbit, y,,(t) is even. Thus equation (11) can be reduced to ur A( to) = -26 y;(t) dt+27 cos wto y,,(t) cos wt, dt I0 I0 =
-2&Y +2-7 cos Wt&.,(W),
where Iher= jy yo( t) cos wt dt. Therefore, only if 7’
Yhet
the heteroclinic =
~~lbwr(W).
(14) bifurcations
will occur if and (15)
THREE-WELL
POTENTIAL
OSCILLATOR
4. COMPUTER
23
CHAOS
STUDIES
This section is purely for the purposes of numerical verification of the results obtained in the previous sections. Thus, we will choose the unstable equilibrium position x0 = 0.5 rather than 0.435 for the experimental system. The damping ratio S is chosen to be 0.1 (S = 0.03 in experiments). 4.1.
SIMULATION
OF
THE
BIFURCATION
CRITERIA
FOR
CHAOS
Our algorithms for computing the bifurcation conditions (13) and (15) basically include two steps. We first numerically solve the differential equation (6) for yO(r) using the fourth order Runge-Kutta technique, and then we numerically integrate the integrals (Y,Iho,,, and ZAP,. Even though the time intervals in those integrals are infinity, only finite time intervals are needed in our computations since yO(t)‘s exponentially converge to the corresponding saddle points, at which y,Jco) = 0. For our three-well potential system, the total length of time in the computations was 10 non-dimensional time units for both the homoclinic and heteroclinic orbits. The initial conditions at p1 and pz will give half of the homoclinic orbit and half of the heteroclinic orbit for forward time integrations. Once yO(I) is found, the integrals (Y, Zhomand Z,,e, are simply a matter of finite summations over the time interval O-10. The critical values of bifurcations are then obtained with the use of the formulas (13) and (15) for the given parameter values. Because the integration processes were involved, the time step 0.01 was used in the Runge-Kutta algorithm for computing yo(r)‘s. In other cases, the time step is taken to be 0.1. Note that the critical values for bifurcations, y,,Yhom and yher, are functions of frequency w. Saving the solutions of yO(t)‘s, one only needs to compute the integrals Zho,,,and I,,,,, and then YI,~,~and yher, for each value of w. For x0 = O-5 and S = 0.1, 50 homoclinic values and 50 heteroclinic values of y were computed over a frequency range of O-6-1.2. The results are shown in Figure 3, where the solid curve is for homoclinic bifurcations, while the dashed curve is for heteroclinic bifurcations. In the next few subsections, we will investigate what happens to the solutions as one crosses these bifurcation curves. 4.2. COMPUTATIONS OF THE STABLE AND UNSTABLE MANIFOLDS The computations of the homoclinic or heteroclinic orbits for the unperturbed system are straightforward. They are obtained by simply numerically integrating (6) with the proper starting points on the orbits. However, when the T-periodic perturbations are
o.ootl 0.6
0.7
0.9
0.9
I.0
1.1
I.2
Frequency
Figure
3. Homoclinic
(---)
and heteroclinic
(- - -) bifurcation
curves in the (w, y) plane;
x0 = 0.5, 6 = 0.1.
24
G.
X.
LI AND
F. C.
MOON
added, the computations of the manifolds have to be performed in the Poincare section. In this case the flow takes place in a three-dimensional phase space (x, y, wt) where the phase of the periodic force is taken modulo 27r/w. When the perturbations are sufficiently small, it can be shown that the Poincare map also has saddle points and stable and unstable manifolds which are close to the corresponding ones for the autonomous system (Y =O). To compute the manifolds in the Poincare section, one first has to find the corresponding saddle points. To be specific, suppose x0 is a saddle point for the following autonomous dynamical system i = f(x),
(16)
where x = (x, y) are the state variables, and f(x) is a vector field where f(xo) = 0. When a T-periodic forcing is added, i = f(x) + &g(t),
(17)
we expect that the saddle point will be slightly shifted for E <<1 in the Poincare section. To find the shifted saddle point x,, one solves the solution to (17) and then determines the PoincarC map by letting t = nT, n = 0, 1,2, . . . . In general the global Poincare map in the plane cannot be determined analytically. However, to find x,, one can linearize equation (17) in the neighbourhood of x0 by letting x = x0+ U: ti = Df(xo)u + &g(t) = Au + &g(t),
(18)
where A = Df(x,) is the Jacobian matrix of f(x) at x0 and is a constant matrix. By doing so, the solution to equation (18) can be obtained as I
to) + u(t) = e A(‘-‘,+L4(
E
I
eActp7)g( 7) dr.
(19)
‘0
One obtains the Poincare map in the neighbourhood u(T+to)=eAru(lo)+E
of x0 by letting t = T+ to
T+‘C, e A’T+‘dg(~) I ‘,I
d7,
and, thus, the perturbed saddle point x, =x,+e(Z-e
AT))’ I,““J eA(T+r,,-~)g(r) d7.
(20)
Note that x, is indeed E-close to the unperturbed one x0: 11x,-x0(/ = O(E). Once the perturbed saddle point x, is determined, one linearizes the vector field f(x) at x, and then determines the eigendirections of IIf( which give the tangents to the corresponding stable (with negative eigenvalues) and unstable (with positive eigenvalues) manifolds W” and W” crossing x,. To compute W”, we choose a small segment along the unstable eigendirection centred at x,. Many points along the segment are chosen to be initial conditions for equation (17), which can numerically integrated. The Poincari map is constructed from the solution of equation (17) by storing points which have time at nT (n=0,1,2,. . .). The resulting points lie on the unstable manifold of the Poincare map of the flow (17). The stable manifold W” can be computed in a similar way, with negative time steps. We now show some stable and unstable manifolds for system (5) in the Poincare map section (see Figures 4 and 5). We want to verify whether or not the stable manifold W” meets the unstable manifold in the Poincari map plane when the driving amplitude reaches the bifurcation value for which A = 0 in equation (12). When x0 =O*S, the
THREE-WELL
Figure 4. Homoclinic
bifurcations
POTENTIAL
OSCILLATOR
of system (4) in a Poincark
section;
CHAOS
x0 = 0.5,s
= 0. I, w = 1.O and y = 0.0334.
t Figure 5. Heteroclinic
bifurcations
of system (4) in a Poincark
section;
x0 = 0.5, 6 = 0.1, w = 10. and 7 = 0.12.
equilibrium point (x,, yO) = (0.5, O-O) is a saddle point for the unforced system (y = 0). For the fixed parameter values 8 = 0.1, w = 1.0, the saddle point for the Poincare map is moved to (0*476,0*00176) at the homoclinic bifurcation value y=O*O334, and to (0.4154,0-00615) at the heteroclinic bifurcation value y =O-117. As shown in Figure 4 for the homoclinic bifurcations, W’ just meets W” tangentially at a few points at the predicted critical value of y = OeO334.For the heteroclinic bifurcations (see Figure 5), it appears that when y = O-12, w” and W” meet, compared with the predicted value of O-117. Note that while the homoclinic bifurcation occurs at y = O-0334 in Figure 4, the stable and unstable manifolds, corresponding to the heteroclinic orbits, are still far away from each other.
G. X.
26 4.3.
LYAPUNOV
LI
AND
F. C.
MOON
EXPONENTS
Having obtained the critical values of the homoclinic and heteroclinic bifurcations, we are in a position to search for the system chaos. We want to investigate, if a homociinic or heteroclinic bifurcation occurs in the system, what will happen with the system response. We answer this question by employing the Lyapunov exponent technique in the driving frequency and amplitude plane (see, e.g. reference [12]). If a motion is chaotic, then all neighboring motions will diverge from it exponentially. The Lyapunov exponent is simply a constant measuring such divergence or convergence on average. For a two-dimensional dynamical system, like system (5), there are two Lypunov exponents associated with each solution. If at least one exponent is positive, then the motion is said to be chaotic. For our purpose of finding chaos, only the largest exponent is needed and simulated for each motion. The advantage of the Lyapunov exponent technique in the study of chaos is that one can systematically investigate the chaotic behavior of the system in a region in a parameter space by simulating the motion and its Lyapunov exponent simultaneously. For the dynamical system (5), we calculated 100 x 100 Lyapunov exponents in the (w, y) plane for the initial condition (0.5,O.l). Only those parameter values of (w, y) with positive exponents were stored and plotted. In Figure 6, each + represents a positive exponent with value greater than, or equal to 0.01 at the given driving frequency and amplitude. The simulation result immediately tells us that many, but not all, parameter values above these bifurcation curves do cause chaotic responses in the system. The higher the driving amplitude is, the more likely chaos is found. We also carried out similar simulations for the initial condition (0.49,0.0), and the result was found to be qualitatively the same as the one in Figure 6.
++
J+
+:+
+
+
+
&,A
i’+*jd+++++ / + + ++ + + +
+ 0.08 -
/
/
+-
/
/
/
/
+,
++/ j4+ t + + + :+$+ / ++
+ t
/’
0.06 -
*++1
//+ ++,#
+
_
+ +
+
/’
+ _ +
_
, /’ 0.04 >./’ _ + 0.02 -
Figure 6. Positive Lyapunov bifurcation curves: 8 = 0.1.
t
.H’ +
exponents
in the (w. y) plane;
chaotic
motions
become
possible
above
the
The Lyapunov exponents in Figure 6 were simulated in a Micro VAX computer. The integration step size was 0.1, and each exponent was averaged in a time period of 650 non-dimensional time units. By checking some of those positive exponents marked by + in Figure 6, we found that they were indeed positive within the time period. However, if the averaging time period is larger than 650, no guarantee can be given to their final values. In fact, experience tells us that an apparent chaotic motion sometimes escapes to a regular pattern after a sufficiently long term. In the next subsection, we present two
THREE-WELL
POTENTIAL
OSCILLATOR
27
CHAOS
phase portraits for some special parameter values in Figure 6 to see the phenomenon transient chaos and persistent chaos. 4.4.
PHASE
of
TRAJECTORIES
To see the effect of the homoclinic and heteroclinic bifurcations, several parameter values of (w, y) were chosen across the bifurcation curves in Figure 6, and the phase portraits were simulated for damping 6 =O*l and the initial condition of (x0, y(,) = (0.5,O.l). For (w, y) = (O-918,0*0435), which is just above the homoclinic bifurcation value of (0.918,0*0312), the phase trajectory is shown in Figure 7(a). The time period for the simulation was 600 non-dimensional time units. We found that the motion appeared random-like in the first 200 time units and then gradually converged to a limit cycle. We also monitored the Lyapunov exponent simulation together with the phase trajectory, and found that the exponent converged approximately to 0.12 for the first 200 time units and then slowly decayed toward zero. In the literature, such a motion is referred to as transient
chaos.
To find a permanent chaos in the time period of 600 time units, we chose another point, (O-714,0*132) from the middle of the cloud in the (w, y) plane (see Figure 6), for simulating another phase protrait. In Figure 7(b) is shown the phase trajectory for this case, and one can clearly see that the motion was indeed chaotic within this time period. Therefore, for our three-well potential system, numerical evidence in Figure 7 suggests that both the homoclinic and heteroclinic bifurcations are necessary conditions for the
r. 5 9
0.4
0.6
0.5
0.7
0.6
-0.6 -1.5
-1.0
-0.5
0.0
0.9
I.0
0.5
1.1
I.0
1.2
I.5
Disdacement
Figure 7. Phase trajectories; (a) transient chaos, 6 = 0.1, w = 0.918, y = 0.0435; (b) persistent chaos, S = 0.1, w =0.714, y =0.132. The simulations were carried out in a period of 600 time units with initial condition at (0.5,O.l) for both cases.
28
G. X. LI AND
F. C. MOON
occurrence of long-term chaotic motions of the system. Satisfaction of conditions for at least one of these bifurcation appears necessary for transient chaos to occur. 4.5.
FRACTAL
BASIN
BOUNDARIES
Basin boundary simulations are another means for checking homoclinic or heteroclinic bifurcations [9, lo]. When a homoclinic bifurcation occurs, the stable manifold W” will (a)
Figure
8. Basin boundaries:
(a) smooth,
y =0.02;
(b) fractal,
y =0.034;
(c) more fractal,
y=O.O5;
S =O.l.
THREE-WELL
POTENTIAL
OSCILLATOR
CHAOS
29
intersect with the unstable manifold W” infinitely many times in the Poincare section. Suppose that the dynamical system has two or more regular attractors, such as two limit cycles for example. A choice of two initial conditions on either side of a stable manifold will lead to two different final states of motions. In fact, the basin boundary in the initial condition plane is the closure of the stable manifolds of the system, just like the Poincare map is the closure of the unstable manifolds. Homoclinic or heteroclinic bifurcations lead to complicated intersection of the stable and unstable manifolds. Therefore, at these bifurcations the basin boundaries will become extremely complicated too, and are called fructal basin boundaries. Such fractal boundaries indicate that whether the system is attracted to one or or the other periodic attractor may be very sensitive to initial conditions. Thus small uncertainties in the initial conditions can lead to unpredictability of the system output even when the motion is not chaotic. This sensitivity to initial conditions, however, may become a precursor to chaotic dynamics. Numerical simulation was carried out in order to determine the basin boundaries for the three-well potential system (5) for different values of driving amplitudes y at S = 0.1, o = 1*O(see Figure 8). The system has three attractors in the neighborhoods of x = 0, * 1. For small y, the steady state motion will be a periodic limit cycle in one of the three attractors. Even if y is increased beyond the critical values of bifurcations, it is still possible that the final steady motion could be periodic rather than chaotic. In fact, near the bifurcation curves (see Figure 6), no permanent chaos was found in our simulations. The colors blue, yellow, and red were used in our simulations to denote initial conditions (starting points), the final steady states of which converged to the left (x = l), middle (x = 0), and right (x = 1) attractors. In Figure 8 are the black and white versions of the original color simulations. In Figure 8(a), y =0.02, which is smaller than the critical value 0.0334 for the homoclinic bifurcations. The boundaries are clearly smooth. Figure 8(b) corresponds to the predicted critical value 0.0334 (= yhom) at w = 1.0. Here many sections of the basin boundaries become non-smooth or fractal. If initial conditions are chosen near the fractal basin boundaries with a small uncertainty in this case, then the final states to which the motion is attracted is not predictable. The uncertainty of final states in the presence of small variations in the initial conditions is a feature of chaotic systems. With further increases in y, the fractal basin boundaries spread out and become more complex, as shown in Figure 8(c) for y = 0.05. By now it is clear that a homoclinic or heteroclinic bifurcation is a necessary condition for chaotic responses of our three-well potential system and a sufficient condition for the appearance of fractal basin boundaries. Numerical simulations tell us that the probability of finding chaos is higher for reasonably larger driving amplitudes.
5. EXPERIMENTAL CHAOS OF THE NONLINEAR SYSTEM In this section we investigate the chaotic responses for the experimental system. Our purposes are twofold: (1) to see if chaos exists in the system; and (2) if it exists, to see if the mathematical model (5) can be used to predict chaos. The external forcing of the experimental system was sinusoidal, supplied by an electric shaker. To measure the motion of the beam, two strain gages were mounted on the upper end of the beam near the clamped end, which produced electrical signals proportional to the bending of the beam. The velocity signals were obtained with the use of an electronic differentiator. The experimental set-up is shown in Figure 9. Without the magnets, the lowest three natural frequencies of the beam were measured to be 4.6, 26.6 and 73.6 Hz. When the magnets were added, the beam had three stable equilibrium positions. The lowest natural frequencies of the beam at the three potential
30
C.
Figure 9. Schematic
X.
LI
AND
F. C. MOON
of the data acquisition
system
of experiments.
wells were found to be 7.8, 8.1 and 8.3 Hz in the left, middle, and right potential wells respectively; they were not symmetric as expected. The asymmetry might be due to prestress in the beam and/or different magnetic strength in the magnets. The experiments were performed at the driving frequencies from 5 to 10 Hz, which were close to the three lowest natural frequencies of the beam. When the driving amplitude was small for a fixed frequency, the response of the beam was periodic of period one. However, when the driving amplitude was increased, we often observed subharmonic motions with period two or more. When the amplitude was increased beyond some critical value, the response would become chaotic. In our experimental studies, Poincare maps were used as a means to distinguish chaotic signals from periodic motions. If a motion was periodic, the Poincare map in the phase plane contained a finite number of points. When the motion was chaotic, there appeared a fractal-looking Poincare map in the (x, li) plane. This highly ordered, fractal-like structure is a basic feature distinguishing chaos from random noise, and is often called a strange attractor. The procedure to experimentally construct a Poincare map can be found in a book by the second author [12]. For the three-well potential system, a typical Poincare map for a chaotic motion is shown in Figure 10. For the purpose of comparison, a numerical Poincare map was also simulated by using equation (3). The parameter values used were /3 =0.577, x0=0.435, w =0.833, which correspond to the experiment, and 6 = 0.15 and f = 0.36, which were higher than the measured ones of S = 0.03 and f= 0.08. We should point out that it was difficult to generate numerically a Poincare map for small damping 6. When 6 is small, it is difficult to find a persistent chaotic response (long-time simulation) which is needed to generate a Poincare map with sufficiently many points. Also, when 6 is small, the Poincare map is very diffuse (close to a forced Hamiltonian system) and hardly any fractal pattern can be observed. In Figure 11, 4000 points (corresponding to 4000 forcing periods) were generated. To find regions in the driving frequency-amplitude plane in which motions of the dynamical system are chaotic, the driving amplitude was slowly increased for a fixed driving frequency near the natural frequencies. When the Poincari map of a response became fractal-looking, like the one in Figure 10, then the motion was declared to be
THREE-WELL
Figure oscillator
10. An experimental is 0.03.
PoincarC
POTENTIAL
map of chaotic
motion;
31
CHAOS
OSCILLATOR
the damping
ratio measured
from the physical
0.6-
2, 0.4.E s 9
0.2-
o.o-
-0.2-
-0.4~""""""""""""""~7 -1.5 -1.0 -0.5
0.0
0.5
I.0
1.5
Dlsplocement Figure 11. A PoincarC w = 0.833, j’= 0.36.
map
from
numerical
simulation;
in equation
(3), 6 =0.15,
/3 =0.577,
x,=0.435,
chaotic at that particular parameter value. Repeating the above procedure for different driving frequencies in the range of 5-10 Hz, we obtained an experimental threshold for chaos in the driving frequency and amplitude plane (see Figure 12). We point out that all motions in the experiments had initial conditions in the left potential well. The duration of observation of a chaotic motion was about four minutes. We also should point out that no guarantee can be given whether or not those chaotic motions are still chaotic at a later time. Also, if a response of the system is chaotic for a particular initial condition, it may not be chaotic for other initial conditions. In Figure 12, we also plot the corresponding homoclinic (dashed curve) and heteroclinic (solid curve) bifurcations by converting the non-dimensional quantities into their physical counterparts. The homoclinic bifurcation curve is considerably lower than the experimentally measured one. The heteroclinic one seems much better, especially at higher frequency. It is clear that while the analytical criteria (13) and (15) are satisfied, these
32
G.
X.
LI
AND
F.
C.
MOON
121,,,,[,,,,,,(,,
,(,,
,,,,
i
\
----o~-~~~.t..~.-c~.5
---__
__
_ -II--%__ 6
7
6
Drlvmg
-----_ I 9
IO
frequency
Figure 12. Criteria for chaos for the three-well potential dynamical formula (22); -, heteroclinic bifurcation (15); - - -, homoclinic
system. bifurcation
-, Experimental;
-
-, heuristic
(13).
lower bounds lie too far below the experimental results. In the following, we derive another criterion from a different point of view. Following an idea introduced in the study of the two-wall potential system by the second author [2], we also calculated a similar heuristic criterion to predict chaos for the three-well potential oscillator. The basic idea is to look for when periodic orbits become large enough to jump out of one of the potential wells (see also reference [12]). As mentioned earlier, the response of the system is periodic for small external forcing. A relation between the driving frequency and amplitude, also involving the response amplitude, could be found through the periodic-solution-solving process for small forcing. Based on experimental observations, we found that when the phase plane orbit become large enough, about the size of the unforced separatrices (the homoclinic and heteroclinic orbits in Figure 2), the motion often had a strange behavior leading to chaos. This property of the solutions was used to obtain a criterion for predicting chaos similar to the homoclinic/heteroclinic criteria (13) and (15). Since all chaotic motions in our experiments had initial conditions in the left potential well, we rewrite equation (3) about the left potential minimum, x = -1. This can be accomplished by replacing x with x - 1, x+&i+
f
;=,
aix’=fw2COS(Wt+#J),
(21)
where (Yr=2P(l
-xi),
a2 = -/?(7 -3x&
a3
=
P (9
-
d,,
a4
=
-sp,
(Yg=p.
Here a phase difference 4 has been introduced for the convenience of mathematical manipulations. In the master’s thesis of the first author [ 111, the harmonic method and the Lindstedt singular perturbation method have been used for the purpose to solve for an approximate periodic solution of the system under the assumption that f is small and w is close to the non-dimensional natural frequency. Doing so, we obtained f=(a/~~){~~~~+[&~~+(~(~~-5c&6cz,)u*]~}”~
(22a)
from the harmonic balance method, and (22b)
THREE-WELL
POTENTIAL
OSCILLATOR
CHAOS
33
from the singular perturbation method. In equations (22), w0 = J2p( 1 -x$) is the normalized natural frequency at the left potential well, and a is the amplitude of the fundamental component of the system response. Notice that both of equations (22) were derived under the assumption that 1w - q,I cc 1, and in equation (22b) we assumed x = a,+ a cos wt + a, cos 2wl, where a, 2 ---~cY~u~/cY,and a, = $x2a2/cx, following the analysis of the relative orders of magnitude. Obviously, under this assumption, the two equations are almost identical. When w is tuned away from wO, the error will be expected. The reason for deriving two relations at the same time is that when w is not close to w,,, probably each of the two relations is good for certain frequency ranges. This is the case in our analysis; the driving frequency can be quite far away from the natural frequencies. This criterion is based on the heuristic idea that chaos might occur when the velocity of the forced periodic motion attains a critical value. If a is a measure of the amplitude of motion, then the critical velocity is chosen to make the orbit in the phase plane close to the size of the separatrix of the unforced, undamped motion. This criterion proved successful in the case of the two-well potential problem (see reference [I]). To represent the closeness of the response to the separatrices, we use the velocity of the beam as a parameter, which can be approximated by aw(=u). If we once again use a “cheat” factor (Y(O< (Y< l), the above idea can be qualified by U=(YU”
or
wa = au0
(23)
near the chaotic regions, where z10is the maximum velocity along the separatrices for the unperturbed system. Substituting it into equations (22), we obtain relations between w and f, with LYas an adjustable constant to match with the experimental data. We found that by properly choosing the factor (~(=0.87), the relation (22a) compares well with the experimental data for w < wO, while the relation (22b) compares well with the data for w > wg. In Figure 12, these two curves are matched at the natural frequency w0 and give a much better prediction to the region of chaotic motions.
6. CONCLUSIONS that homoclinic or The results in this paper, once again, lead to the conclusion heteroclinic bifurcations in a dynamical system can be accurately predicted with the use of Melnikov’s method (see Figures 4 and 5). Chaotic motions of the three-well potential system were found both experimentally and numerically. A horseshoe map is found for the system when homoclinic or heteroclinic bifurcations occur, which often causes complex dynamics. It is important to be aware that although the existence of horseshoes gives the possibility of chaos, one cannot conclude the existence of chaos without additional information. The Poincare map technique is the best means to check the existence of a strange attractor motion in a laboratory. If the motion is really chaotic, then a highly ordered, fractal-looking, pattern should be observed. On the other hand, the Lyapunov exponent technique is very useful in checking chaos in numerical simulations. When doing so, care must be taken regarding the convergence of the exponent. The three-well potential system has three limit cycle attractors for sufficiently small forcing. The existence of homoclinic or heteroclinic bifurcations implies the existence of fractal basin boundaries of attracting domains and the loss of absolute predictability in the face of small uncertainties in initial conditions. In fact, the Melnikov technique is accurate in predicting the occurrence of such fractal basin boundaries. The experiments of the three-well potential system also showed chaotic motions, even when the forcing inputs were deterministic. In the frequency-amplitude plane, the experimental criterion for chaos is above both the homoclinic and heteroclinic bifurcation
G. X. LI
34
AND
F. C. MOON
curves, indicating that these bifurcations are necessary conditions for the appearance of chaos in the system. The heuristic criterion (22), however, gives a much closer prediction
for chaos.
REFERENCES
of Applied Mechanics 47, 638-644. Experiments on chaotic motions of a forced nonlinear oscillator: strange attractor. 2. P. J. HOLMES 1979 Philosophical Transactions of the Royal Society A (London) 292, 419-448. 1. F. C. MOON 1980 Journa/
A nonlinear oscillator with a strange attractor. 3. V. K. MELNIKOV 1963 Transacrions of the Moscow Mathematical
4. 5. 6. 7. 8. 9. 10.
Sociefy 12, l-57. On the stability of the center for time periodic motion. F. C. SALAM 1986 SIAM Journal of Applied Mathematics 47 (2), 232-243. The Melnikov technique for highly dissipative systems. F. C. MOON, J. CUSUMANO and P. J. HOLMES 1987 Physica 24D, 383-390. Evidence for homoclinic orbits as a precursor to chaos in a magnetic pendulum. J. GUCKENHEIMER and P. HOLMES 1983 Nonlinear Oscillafions, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag. A. WOLF, J. Swrrr, H. SWINNEY and J. VASTANO 1985 Physica 16D, 285-317. Determining Lyapunov exponents from a time series. A. J. LICHTENBERG and M. A. LIHERMAN 1983 Regular and Stochastic Motion. New York: Springer-Verlag. F. C. MOON and G. X. LI 1985 Physica/ Review Letters 55, 1439-1442. Fractal basin boundaries and homoclinic orbits for periodic motion in a two-well potential. C. GREBOGI, S. MACDONALI), E. OTT and J. YORKE 1983 Phvsics Letters 99A, 415-418.
Final state sensitivity: an observation to predictability. 11. G. X. LI 1984 M. S. Thesis, Cornell University, Ithaca, N. Y. Chaotic System with Five Equilibrium States. 12. F. C. MOON 1987 Chaotic Vibrations. New York: John Wiley.
Vibrations
of a Nonlinear