Study of the Chaotic Behavior of a Damped and Driven Oscillator in an

Study of the Chaotic Behavior of a Damped and Driven Oscillator in an Asymmetric Double-Well Potential Sandra L. ~ a e e , Carole ' L. LeMaster, and Clifford B. ~ e ~ a s t e ? Boise State University, Boise, ID 83725 Donald A. McQuarrie University of California, Davis, CA 95616 Very few real systems demonstrate regular predictable behavior and the mnionty of thode that do are determined by an initial set of cmdkions that produce a single solution. Many nonlinear dissipative systems exhibit behavior that is deterministic, but unpredictable through time. The study of these systems is now called the study of "chaos"or "nonlinear dynamics". This new field has been brought to the attention of the general public by the best-selling book entitled Chaos, Making a New Science by James Gleick as well as numerous other books (14). Despite the profound implications of chaos and its widespread appearance in many fields, including biology, chemistry, economics, engineering, medicine, meteorology, and physics, undergraduate students are seldom introduced to its existence, let alone its study. This leads them to believe that all systems are modeled by equations having closed form solutions and that their behavior can be ~redictedfor long periods of time. Data not supporting these conclusions are rejected as aberrant. Because of the current accessibility of reasonably powerful computers, it is no longer necessary for students to confine their studies to systems having only analytical solutions. Most science and engineering majors, who have reached their junior year, have the mathematics, physics, and c o m p u t e ~ p r o g r a m ~ skills n g necessary to tackle chaotic problems. A wealth of chaotic systems suitable for study have been described (5-7). Projects can be scaled to fit any time frame by limiting either the number of parameters to be varied or the number of methods used to model behavior. The study we describe here was completed in two semesters as a senior undergraduate research project. Chaotic behavior results when inaccuracies or uncertainties in the starting conditions grow exponentially, as oooosed to linearlv. with time: conseauentlv. the svstem's state shortly becomes totally &predic&ble."This condition reflects a n extraordinarilv sensitive de~endenceon the oarameters defining the sykem's initialitate. Consider the following quote from an article by Grebogi, Ott, and Yorke

(8). According to Laplaee, determination of the future depends only on the present state. Chaos adds a basic new aspect to this mle: small errors in our knowledge can grow exponentially with time, thus making the long-term prediction of the future impossible. All that is required for a system to be potentially chaotic are one or more nonlinear coupling terms in its equations of motion and more than two phase-space dimensions. Chaotic systems are generally modeled by a set of three or A Drelirninaw version of this article was Dresented at the 203rd Na-

more nonlinear differential equations. The damped and driven oscillator in a n asymmetric double-well potential is such a system, having three phase-space dimensions (position, velocity, and time) corresponding to one and a half degrees of freedom and a nonlinear (cosine) term in its equations of motion, thus satisfying the conditions necessary to produce chaotic behavior. Asubset of this system, a damped and driven oscillator in a symmetric double-well potential, the "Duffing Oscillator", has been discussed previously (7).We proposed to study the effect of asymmetry on this system's behavior. An oscillating particle in a double-well potential can serve as a model for several physical and chemical systems. Holmes (9)first used the DuiXng equation to model the forced vibrations of a buckled beam. An elastic magnetic beam is suspended from a cantilever support. Two magnets are placed near the free end of the beam, creating two stable equilibrium positions. The entire apparatus is then placed on a vibrating table that oscillates with a certain amplitude and frequency. For small amplitudes, the beam vibrates about one of the equilibrium positions, but as the amplitude is increased, the beam can jump out of the potential well and move chaotically from one well to the other. If the magnets are of different strengths, the potential function is asymmetric. Several chemical systems that can be described by a n oscillator in a double-well potential become apparent if we understand the physical significance of the model a s illustrated by the buckled beam. Simply, it describes the motion of a particle forced between two stable equilibrium positions. If the potential energies of the positions are not equal, then the potential function is asymmetric. The classical motion of a hydrogen atom moving between two potential hvdrozen-bond in^ sites within a molecule (as in polypeptides ind proteins) or between two sites m ifferent molecules (as in the case of DNA, can be described in such a way. he rotational motion of groups of atoms, about molecular bonds, between two stable conformations (i.e., the rotation about the Cz-Cs bond in n-butane) is described by torsional potential energy functions that resemble the double-well potential. The anti andgauche forms of n-butane lie at the bottom of energy wells separated by -4 kJ mol-'; thus, the potential energy function is asymmetric. Diffusion through doped (impure) solids and oscillations in plasmas are two other examples. Avery active area for the study of chaos in chemistry is reaction kinetics. Any reaction process described by a system of three or more kinetic rate equations coupled by terms depending on the concentration of more than one reactant is potentially chaotic. An excellent description of the theoretical and experimental techniques used to study chaotic behavior in eas- and solution-ohase reactions. heterogeneous catalysis, electrodissolutiou reactions, and biochemical svstems is eiven bv Steohen Scott in his monograph entiiled ~hemical~hrEos(lb).

-

'Author to whom correspondence should be addressed

Volume 71 Number 5 May 1994

363

Chaotic behavior in the biochemical processes described above: intermolecular hydrogen-bonding and peptide-bond rotation in proteins, intramolecular hydrogen bonding in DNA, and in biochemical reactions may relate to the causes of certain diseases and aging. An understanding of the conditions that bring about chaotic behavior may eventually lead to treatments for and prevention of these conditions.

U"(l)=!k+6d+12e+20f+30g>0 (2) W(-l)=2~-6d+12e-20f+30g>O (3) Substitution of the solutions for the constants a through f into the inequalities (eqs 1,2,and 3) and solving the inequalities forg yields the following conditions forg.

Definition of the System

Thus,

The system we describe is that of a particle, whose motion is conhed by a double-well potential, driven by a periodic external force of constant amplitude and subject to a frictional drag proportional to its velocity. Below, we discuss, in turn, the various forces acting on the particle. The Asymmetric Double- Well Potential

The asymmetric double-well potential, Ux), can be described by a sixth degree polynomial. ~ ( x ) = +bx+a2+dz3+ex4+fx5+ffr6 n The with the constants, a through f: 'peeified as two wells of the potential are chosen to lie at x = f1 and the intervening maximum at x = 0. Therefore, we write atx= 1, U(1)=0 righthand well ,eR-hand well at x = -1, U(-1) = $ at x = 0, U(0)= a maximum The parameters P and a correspond to the heights at x = -1 and x = 0, respectively. These values of x describe equilibrium points (minima and maxima) on the potential; therefore, the first derivative of the potential function at these points is zero and we write atx=l, LP(l)=O right-hand well atx = -1, U'-1) = 0 left-handwell at x = 0, WO) = 0 maximum This yields the following set of six equations: U(0) = a = a U'(0) = b = 0 U(l)=a+b+e+d+e+f+g=O V(l)=b+!k+3d+k+5f+6g=O U(-D=a-b+c-d+e-f +g=P U'(-l)=b-!k+3d-k+5f -6g=O

g>$-a

g>%-a

U(x)=a+(p-2a+g)r with ~ a - p > g 16 >g-a Each of the three parameters, a,P, andg, have physical ,ignifiCance directly related to the graphical representstion of the potential. The height of the unstable equilibrium position at x = 0 is a.The asymmetry parameter, P, is the height of the left-hand well. Finally, the value ofg corresponds to the broadness of the saddle at the unstable equilibrium position and the widths of the two wells, Recall thatg is restricted by the values of the second derivatives that define the direction of curvature of the equilibrium positions. Figure 1illustrates the effect of variations in p andg on the shape of the double-well potential. The Equation of Motion for the Oscillator

The equation of motion for the oscillator is derived from the potential, a damping term proportional to the velocity, and a harmonic forcing term. -d2x

dU

dt2 -

dx potmtial

dx -e +Few (wt) dt

damping

(4)

fomieing

where U is the potential, x is the position, t is the time, c is the damping constant, F is the forcing amplitude, w is the forcing frequency, and, finally, fix)=--dU dx 023

which, when solved simultaneously for t h e constants, gives a=a, b=O a a c=P-Za+g d =9 4 -$+2a-* 32 e= f= 2

-

,

gxand k e e ~ i n the e remainder). To keep Afa small enoueh to ensure acc;ra& of the solutions, but iarge enough sothat excessive computer time is not used, it is adjusted after each calculation based on the difference, h z ~If. the difference is greater than a predetermined value, E ,Afa is increased and, if less than E, it is decreased for the next calculation. We used the increase and decrease amounts suggested by and 0.95 (dhz2)1'5,reBaker and Gollub (12),0.95 spectively. - o u r trajectories were generated a t specific values of the parameters c, F, and w. Typically, transient hehavior was A.

4.6

Position Figure 2. Two orbits of the damped and driven oscillator in an asymmetric double well potential calculated using standard parameters. but started from slightly different initial conditions. The trajectory represented by the solid line was started at xl = 4 . 1 , x, = 0.5,x3=0 and by the dashed line at xl = 0, x, = 0.5, x3 = 0. The orbits diverge exponentially, characteristic of the sensitive dependence on initial conditions demonstrated by chaotic systems.


365

-1

99

199 299 399 Number of Drive Periods

499

Figure 3. The convergence of the three Lyapunov exponents for the oscillator calculated fora double-well potential with a = 0.25. I3 = 0.1, and g = 0, a damping constant of 0.5, a forcing amplitude of 0.5. and a forcingfrequency of 1. Note the positive exponent indicating a region ofchaotic behavior. A positive exponent of 1 bit per drive period. as shown, means that the solution calculated from initial conditions accurate to 10 ppm (-23 bits) will become unpredictable, except to say that it lies on the attractor, in 2311 or 23 drive periods. removed from the system by ignoring the solutions for the first 12 drive cycles (75time units, 2n time units per drive ueriod). This value was increased as necessaw for some brbits.The time stepsize was initially set at 0.2513 (25 increments per drive cvcle) and adiusted as described above. ~rajectdries(orbits) may be characterized by their periods. A period-1 orbit closes after one drive cycle, a period-2 orbit after two drive cycles, etc. A chaotic orbit does not close. Poincare Sections Poincad sections are plots of position versus velocity at a specific time or phase anple, 0,in the drive cycle.3 They ha;@ one less dimension than the corresponding phase diakwam; thus, they are utilized to simplify phase diagvams when they b&o&e complicated, as in resons of chaotic motion. For example, a periodic orbit, whose phase diagram is a limit cycle (one dimensional) has a Poincarb Section that is a point (zero-dimensional).The period of the orbit can be determined from the number of points on the section, one point for period-1, two for period-2, etc. Poincarb Sections of chaotic orbits appear layered, have fractal (noninteger) dimension, and are referred to as "strange attractors". They consist of an infinite number of points, which, when repeatedly magnified, reveal an infinite number of layers resembling the overall structure of the attractor. This property of strange attractors is called selfsimilarity. Poincar6 sections are computed similarly to trajectories, except that, each time the solution is calculated, the time, x3 (modulo 2n), for that solution is compared to a specified moment (phase angle) in the drive cycle. If it equals this value, it is plotted on the screen and recorded in a file. In reality, it would be difficult to set the time stepsize accurately enough to repeatedly duplicate this exact moment for every drive cycle. For this reason, successive values of xa are compared to see if the precise moment lies on the intewal between them. If so, time is set to that moment and the solution recalculated and plotted. %me modulo 2n is sometimes called the phase angle. 2n time unils are required to complete one drive cycle: thus, at 3n time units. the phase angle is n or 180. 366

Journal of Chemical Education

1 " "

-1.5

I

i " "

-1

" "

4.5

I " "

0

1 " "

0.5

" "

1

1.5,

Position Figure 4. Concentric orbits generated when the force is varied from 0.01011 to 0.34011 in steps of 0.03. Standard parameters and initial conditions IC1 and IC2 were used to produce the orbits. Our Poincarb sections consist ofjust over 2200 points abstracted (at the specified phase angle) from sets of about 30,000 points. The first 16 drive cycles (100 time units) were ignored as transient and the time stepsize was 0.5026. Bifurcation Diagrams While trajectory plots and Poincarb sections allow us to look at the behavior of the oscillator at specific values of the parameters, bifurcation diagrams provide information about the evolution from periodic to chaotic behavior. A bifurcation occurs when the number of solutions to the equations of motion change as some dynamical variable, such as the amplitude of the driving force or the damping constant is varied. Bifurcation diagrams are generated similarly to Poincar6 sections, except that the whole procedure used for calculating Poincarb sections is placed within a loop, within which the dynamical variable to be varied is incremented from a specified initial to a specified final value. While Poincar6 sections are plots of velocity versus position at a fixed phase of the drive cycle for fured values of the dynamical variables, bifurcation diagrams are plots of either position or velocity versus a dynamical variable. Like for Poincarb sections, the velocity must be calculated at a consistent phase of the drive cycle. For most of the bifurcation diagrams shown below velocity is plotted versus the forcing amplitude. For each value o f F the system was allowed to reach a steady state by ignoring the solutions for the first 40 drive cycles. The velocity at the start (@= 0) of the next 40 drive cycles was plotted. This generates 40 values of the velocity, which are all identical for a period-1 orbit, yield two separate values for a period-2 orbit, etc. The integration time increment used was 0.2513; therefore, each drive cycle (2n)required 25 integrations of the motion equation. A typical bifwcation diagram required 8-10 hours of CPU time on an IBM-compatible 33-MHz 466 personal computer.

Figure 5. Orbits generated from three different sets of initial conditions. The orbit in the left-hand well is from ICI, the outer right-hand orbit from IC2, and the inner right-hand orbit from IC3. Standard parameters were used in the equations. (a) F = 0.3450 and (b) F = 0.3520.

b

Position

F gure 6. (a) Three-dimensional chaotlc traleclory and 1s pro,ecton onto the pos~tion-veloctyplane at F = 0.365. (b) Pomcarb sectlon of ttis tra.ectory. Slanoaro parameters were ~ s e dn the equal ons. lndial conditions are IC2.

Lyapunov Exponents

Lva~unovexwnents are calculated as the mowth rate with &meofthe volume ofan infinitesimal sph;re ofinitial conditions with the initial point on the attractor at thecenter (13). There is an exponent, h, corresponding to each principal axis of the sphere. The exponents are defined in terms of the length,^, of that axis and ordered from smallest to largest.

The sphere defines a space tangent to the phase space. The action of the flow (defined by the equations of motion) causes this sphere to deform into an ellipsoid. The principal axes of this ellipsoid are determined by the evolution, through the linearized equations of motion, of n orthonorma1 vectors tied to the central trajectory calculated from the nonlinear equations. Numerically this is accomplished by simultaneously integrating the nonlinear equations for some initial condition and the linear equations for n different initial conditions, that define the frame of n orthonorma1 axes vectors. For our system, this requires the

simultaneous integration of 12 equations, three nonlinear. and nine linear. This procedure is repeated for each point as the central trajectory evolves in time. The infinitesimal sphere deforms as each initial vector changes magnitude and direction according to the linearized equations. There are two problems that occur when these 12 equations are numerically integrated. First, if the magnitudes of the linearized vectors are allowed to grow too large or too small, the computer can not accurately represent the numbers. Second, all the linearized vectors tend to move towards the direction of most rapid growth and their orientations appear to be identical (at computer precision). These two problems are prevented by repeated GramSchmidt reorthonormalization (GSR) once each drive oeriod. The GSR procedure is explained in detail by ~ o l f e t al. (131and our BASIC oromam to calculate Lvaounov exN pres&tkd in this ponents is based on ~ ~ ~ F & T R Acode article. The Lyapunov exponents are obtained by summing the logarithm of the magnitudes of the three principal axes vectors calculated each time increment and dividing by the total time. In time, the exponents converge to their true values. Volume 71 Number 5 May 1994

367

Figure 7. At values of the forcing amplitude above -0.39, non-transient periodic orbits no longer depend on inaial condiiions. The first bifurcation cascade occurring after this point begins at F = 0.425. (a)Period-I orbit at F = 0.42. (bj Period-2 orbit at F = 0.43. (Cj Period-4 orbit at F = 0.435. (d) Chaotic trajectory at F = 0.45002. The volume of the ellipsoid grows a t the rate of 2@1+ utwhere hi, the ith Lyapunov exponent, is in bits per second or in bits per drive period if multiplied by 2n. A typical Lyapunov exponent spectrum converges in less than 500 orbits and takes about 15 min using a time step of 0.2513 (25 increments per drive period) on a 33-MHz 486 IBM compatible personal computer. The pattern of convergence is shown in Figure 3. Our system has three Lyapunov exponents, corresponding to each of its three phase-space dimensions.Their values provide information about the system's behavior. One exponent, that corresponding to the time dimension will always be zero, because its axis runs parallel to the flow and does not contribute to expansion or contraction of the ellipsoid volume. Contracting axes have negative exponents and expanding axes have positive exponents; thus, a Lyapunov exponent greater than zero is representative of chaotic behavior. Two equal and negative exponents indiorh;t and two aero exponents indicate cate a stablc a chanec in the number of solutions; thus, the occurrence of a bif&ation. +

368


Lyapunov exponents may be used to determine the fractal dimension of strange attractors using eq 8, the KaplanYorke relation (141. h, + h, +...+ hj fractal dimension =j +

1

hj+.jl

(8)

where the exponents, k, are ordered from largest to smallest and j is the index of the smallest exponent that is 2 0. Our reported Lyapunov exponents are the average of the exponents calculated for the 500-599th drive cycles, 25 exponents per drive cycle (At = 0.2513). Reported uncertainties are at the -95% confidence level (20). The reported errors in the fractal dimensions are propagated from the uncertainties in the exponents. Linearization of the Equations of Motion

Lyapunov exponents are evaluated as the rate of the divergence in phase space of two trajectories beginning at initial conditions that are infmitesimally separated. Let % = (x,, X2, ..., x")

the rate of separation for the nth component is

11

Specifically, this is equivalent to the linear terms of a multivariant Taylor series expansion on each component, XI,xz, ...,x,, which can be written in the following fonn.

afi

afl af, -

ax, a*

dtI")=:-

a

-1.5

ax,

af2 afz

afi

SOC,

ax, a%

..' ax,

171

k"

-..

-

-'

-

.

-1

4.5

0

0.5

1

1.5

~~lltion Figure 8. Poimar6 section of the chaotic orbit at F = 0.45002 (standard parameters, @ = 0).The fractal dimension of this strange atlractor isl.13 (0.03).

&"

.

with the partial derivatives evaluated at Xo= (XI,xz, ...,x,). For our system, defined by eqs 5 , 6 , and 7, the linearized equations are

and

& + u = (xl + €a1, x2+ spra, ..,

X"

+ &")

be vectors describing two sets of initial wnditions separated by the infinitesimal difference w. Each set of initial conditions generates a set of solution curves (trajectories) defined by the system of coupled nonlinear first-order differential equations. For a generalized system

d(%, --

[4a- 28 - 2g + dt

9,(12a-

68 - 24g).:

I

- 1 5 ~ -: 30& &, - c&r2- F sin (x3)k8

Results

-0.2

1

0.1

I 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Driving Force Amplitude Figure 9. Bifurcation diagram of velocity versus driving force amplitude calculated wilh standard equation parameters. The driving force amplitude was incremented by 0.001 between 0.200 and 0.865. Black dots indicate solutions for ICl, red for IC2, and green for IC3. Blue dots represent solutions in the region where periodic solutions do not depend on initial wndiiions. The velocity values denoted in red have been displaced by -0.125 in the region F = 0.200 - 0.340 for visibility. The true velocity values are identical to the green values.

When the potential parameters are set at a = 0.25, P = 0, andg = 0, the equations ofmotion, eqs 5,6, and 7, reduce to those of the Duffingoscillator. Because of a previous study (7)of the Duff~ngoscillator, which is representative of an oscillator in a symmetric double-well potential, and our desire to study the effect of asymmetryon this system, we chose to use those values for a and g. To introduce asymmetry to the system p was set at 0.10. The oscillator parameters used in reference 7, c = 0.50 and w = 1.0 (one push per drive period), were also used in most of our study. The driving force amplitude, F, was varied. These values of n.. O.,-. e. c, and w will be referred to as standard oarameters. Initially, we looked a t changes in the trajectories with variations in driving force amplitude using the standard parameters and two sets of initial wnditions, one placing the oscillator in the


.

369

Figure 10. (a)The chaotic orbit at F= 0.6600. (b) Period-12 orbii at F=0.6705. (c)Period-6 orbit at F= 0.6715. (d) Period-3orbit at F = 0.6750. Standard parameters were used. leRhand well, xl = -1, x2 = 0, x3 = 0 (called IC1) and one placing it in the right-hand well, XI = 0.5, xz = 0, x3 = 0 (IC2).For a forcing amplitude of zero, the oscillator comes to rest in the bottom of the well. Its attractor is a point. Figure 4 shows the changes in the solutions as the force is incremented from 0.1011 to 0.34011 in steps of 0.03. The projection of the trajectories onto the position-velocity plane (as shown) is a period-1 limit cycle. The oscillations about the equilibrium points increase in amplitude until, a t F = 0.34011, the limit cycle in the righbhand well shows an abrupt increase. This abrupt change in amplitude is called an "amplitude jump" and is characteristic of nonlinear oscillators. If the initial position, XI, of the oscillator is placed a t 1.0 (IC3), instead of 0.5, this jump does not occur. The amplitude of the oscillations generated from IC3 continues to increase smoothly until F 0.380. Between F = 0.34011 and F 0.371, there are two attractors in the right-hand well and one in the left-hand well. The related Duffing oscillator produces four attractors, two in each well, similar to those in our right-hand well and mirrored in the left-hand well. Note (Fig. 1)that changes in P alter

-

370

-


the shape of both wells; therefore, the values of F, which produce similar orbits, are slightly different for the Duffing oscillator than for the asymmetric case (see ref. 7). At F = 0.34100, the first bifurcation occurs. The orbit generated from IC2 becomes period-2, closing after two drive cycles. The bifurcations continue, period-4 at F = 0.34984, period-8 a t F = 0.35226, period-16 a t F = 0.35280, etc., until the system started a t IC2 reaches chaos a t F = 0.35295. This chaotic window persists until F = 0.372, when the two attractors (from IC2 and IC3) in the righthand well converge to one period-1 attractor. Between F = 0.372 and F = 0.390, there are two attractors, one in each well. Finally, a t F > 0.390, all three sets ofinitial conditions generate thr same whits, initially period-1 in the left well. Figure 5 shows thc thrcc attractors at F = 0.3450 and I.'= 0.3320, and Figure 6 shows the chaotic trajectory from IC2 at F = 0.365 in thrcc dimensions alone with its msition-velocity projection and its Poincar6 section (strange attractor). Note that the forcing amplitude is not great enough to push the oscillator over the maximum to the other well; thus, its strange attracto'r is only a subset of that gener-

-

1

1. and the fractal dimension of the strange attractor depicted by the Poincare section i s 1.15 (0.01). As stated above, a t values of F greater than -0.39, the initial conditions no longer affect periodic solutions to the differential equation^.^ The forcing amplitude is now great enough to push the oscillator over the maximum from the right to the left-hand well and, eventually, it visits both wells. Figure 7 shows the 3D trajectories and 2-D projections for the first bifurcation cascade involving both wells, which begins with a bifurcation to periodd a t F = 0.425 and enters a chaotic window a t F = 0.45002. The Poincare section of the stranee attractor for F = 0.45002 is shown in Figure 8. The calcuJ lated Lyapunov exponents of this 1.5 chaotic orbit are 0.649 (0.017),0, and -5.182 (0.017), giving the Position Poincare section a fractal dimenFlgure 11 Pomcarb sect ons forthe per 00-3,per od-6, and penod-12 orblts are shown a ong wltn that sion of 1.13 (0.03). Another method of viewing the of the chaotlc traiectoty for F = 0 6600 Lyap~novexponents calc~latedfor tne strange anractol' are behavior of the system is by cre0 666 (0 013) an0 4 221 (0 01 3) g vfngit a fracta dmens on of 1 13 (0 02) atine a bifurcation diaeram. Abifurcation diagram showing the values of the velocity calculated for driving force amplitudes between 0.200 and 0.865 usine standard parameters is shown Fieure 9. This diaeram is a composite (overlay) orthree bifurcation diagrams, each using a different set of initial conditions, IC1, IC2, and IC3. Its generation took about 30 hours of computer time. Dots in the colors black, red, and green correspond to velocity solutions calculated starting a t IC1, IC2, and IC3, respectively. Blue dots represent velocity values calculated beyond the point where periodic solutions are dependent on the set of initial conditions chosen. The chaotic windows and bifurcations to period2 described above are readily seen on this diagram. Also apparent is a large period-3 window from F = 0.661 to F = 0.756. In this window, perioddoubling occurs as the forcing Driving Force Amplitude amplitude is decreased, rather Figure 12. Bifurcation diagram showing detail In the region of the first bifurcation cascade generated from iC2. Standard parameters were used and the driving forceamplitude was varied fornl 0.3390 to 4 ~ h statement i ~ may create some 0.3600 in steps of 0.0001. confusion, as chaos has been defined as sensitively depending on the innial ated by higher forcing amplitudes. This orbit is still chaconditions. Periodic orbits, however, are not chaotic. Though caicuiated from a potentially chaotic system of equations, the palticuiar otic, having Lyapunov exponents of 0.831 (0.0091, 0, and parameters used in the equations may not produce chaotic behavior. -5.362 (0.009) bits per drive period. The fractal dimension Once transients are removed from the solution cuwe, non-chaotic of its strange attractor, calculated fmm the Lyapunov exsolutions result. There may be many attractors for a specificset of using eq 8 with j = 2 is 2.15 (0.01). the two-¶meters. This set may be composed of attractors of varying, inmensional position-velocityspace used to plot the Poinear6 cludingfractal,dimensions, Regions containing conditionsthat section, the exponent corresponding to time, 0, is lost,j = to identicalattractorsare calledbasins of attraction.

-

-


371

than increased as for the two windows discussed above. Chaotic, period-12, period-6, and period-3 orbits in this region are shown in Figure 10 and the Poinear6 section for the chaotic orbit with those of the periodic orbits overlaid is shown in Figure 11.For this type of trajectory plot (Figures 6,7, and 10)the period can be determined by counting the number of 3-D solution curves emanating from points on the 2-D projections. This diagram (Fig. 9) shows several other areas of sufficient interest for further investigation including a period-5 window around F = 0.64. Bifurcation diagrams are valuable for the determination of accurate bifurcation points. Although bifurcation points can be determined by varying the parameters used for phase diagrams, this soon becomes tedious, as bifurcation points for higher period orbits become .-loser and closer together and re- Figure 13. Biurcation diagram showing the effectsof asymmetry changes on the solutions to the equaquire resolution to five or six dec- tionsof motion. Pis incremented in steps of 0.001 from -0.200 to +0.125. Standard parameters with F= 0.365 and IC2 were used in the calculation. imal places. The bifurcation nointi reoorted above were de. Substitution into eq 9 gives a value of 4.6 (0.2) for the iermined'from high-resolution bifurcation diagrams created for reeions in which ocriod-doubline is occurrine. FieFeigenbaum number, again supporting the period-dou-bling route to chaos for this system. Bifurcation diagrams w e 12 is &ch a diagram: Resolution is%.0001~ in k e ; eion between F = 0.3390 and F = 0.3600. This diaeram of velocity versus driving force, analogous to Figure 9, for :bows the region of the first hifurcation cascade generated p = 0 (symmetriccase) and p = -0.1are shown in Figure 14. from 1C2. The amolitude "iumo and first three bifurcations Asymmetry can be seen to affect the positions and sizes of are clearly seen. the chaotic and periodic windows. The value of P also afThe oeriod-doublim mute to chaos is known to follow fects the amplitude jump, which gets larger as P increases. universal scaling laws. For example, the ratio between No amplitude jump is seen a t P = -0.10; a t P = 0.05, the successive values of F a t the bifurcation points approaches jump covers about 0.1 velocity units increasing to 0.14, a constant called the Feigenbaum number after its discov0.20, and 0.22 velocity units at P = 0, +0.050, and +0.10. erer. This ratio is defmed by the following: Table 1. Verification of the Lyapunov Sum Rule for the Fk - Fk-i - 4.669 201 609 102 9909... lim Damped and Driven Oscillator in an Asymmetric Douk+ Fkii- Fk (9) ble-Well Potential. hiare in bits s"x In 2 where k is the order number of successive bifurcation c F hl ha points. Bifurcation points accurate to 0.00001 were located h" Ch by incrementing F by 0.00001 in each area close to a bifurcation. For our system, using standard parameters, Fl= 0.200 0.370 0.150 4.350 0 -0.2W 0.34100, Fz = 0.34984, F3 = 0.35226, and Fg= 0.35280. Sub0.300 0.470 0.128 4 . 4 2 8 0 -0.300 sequent bifurcation points become difficult to resolve with accuracy; however, even after only four bifurcations, the 0.400 0.470 0.090 4.490 0 -0.400 value calculated for the Feigenbaum number using eq 9, 0.500 0.480 0.087 4.587 0 -0.500 4.48(0.17), is close enough to postulate that our system 0.600 0.540 0.090 4.690 0 -0.600 does, indeed, follow the period-doubling route to chaotic behavior. Bifurcation diagrams are used to map changes in Table 2. Fractal Dimension of Strange Attractors Versus Damping the solutions to the equations of motion as some paConstant rameter is varied. The diagrams described above show the solution changes with variation in forcing c F hr ha Fractal Dimension amplitude. Because of our interest in studying the effect of asymmetry of the double-well on the 4.462 (0.003) 1.35 (0.03) oscillator's behavior, a bifurcation diagram in which 0.30 0.60 0.161 (0.003) 0 is varied was constructed and is shown in Figure OAO OfiO 0.1 02 (0.002) 4.502 (0.002) 1.20 (0.02) 13. Standard parameters and F = 0.365 were used 0.50 0.60 0.086 (0.001) 4.586 (0.001) 1.15 (0.02) in its calculation. One bifurcation cascade is obvious 0.60 0.m not chaotic in this diagram. The bifurcation points are Pi = - 0.70 0.60 0.110 (0.002) 4.810 (0.002) 1.14 (0.02) 0.1579, Pz = -0.0918, p g =-0.0659, and pq= -0.0617. A

372


,

1.3

I

.s 1.0

0.9

.

0.8 -

-

0.7

6

8 -

-

.

0.6 -

-

0.5 -

0.4 -

.

,-

; /

.,:-. .-, ,.;;

0.2 0.1 0.0

-

*;-

0.3 -

.

>I

-

-

_

.

-0.1 -0.2 0.20

-

0.25

0.30

0.35

0.40

0.45

0.50 0.55

,

0.60

0.65

0.70 0.75 0.80 0.85

Driving Force Amplitude

0.30 and F = 0.39. Exponents were calculated from all three sets of initial conditions IC1, IC2. and IC3. IC1 and IC3 yielded identical exponents in this region, with stable period-l orbits a t F < 0.313. Exponents generated from IC2 show F values correspondingto most stable periodic orbits (two equal or nearly equal negative exponents) a t F < 0.331 (period-1) and F = 0.346 (period-2). Bifurcation points (two zero exponents) are apparent a t F = 0.341 (to period2), F = 0.350 (to period-41, and F = 0.353 (to period*) and two chaotic windows, characterized by one positive exponent are seen between F = 0.353 and 0.369 and from F = 0.381 to 0.386. The sum of the Lyapunov exponents is equal to the time-averagedex~onentialdivergenceor convergence of the flow in phasespace.

where the Lyapunov exponents, hi a r e expressed base e (Lyapunov exponents in bits per drive period must be multiplied by in t2112x). Any dissipative dynamical system will have a t least one negative Lyapunov exponent and the sum of the exponents must be negative for the system to be dissipative. The left-hand side of eq 10, called the Lie derivative, equals -c for our system; thus, the sum of the exponents (righbhand side of eq 10) should equal -c. Table 1 verifies this Lyapunov sum rule for our system. Standard parameters, along with the values of c and F given in the table were used to calcuI ' late thevalue exponents. 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 of the damping eonDriving Force Amplitude stant also affects the fractal dimension. When damping i s Figure 14. Bifurcation diagrams using standard Parameten and IC2 for (a) B = -0.10 and (b) P = 0 large, the fractal dimension is (symmetric). small and points on the strange attractor a r e close together. While bifurcation diagrams are plots of some system There is an overall increase in dimension with a decrease variable, i.e., velocity or position, versus some system pain damping (see Table 2). The dimension change appears rameter, Lyapunov exponent spectra are plots of the expoas an increase in the soacine of lavers in the ~oinca>i ~ - set--~ nents versis some system parameter. ~or-thisreason, they tion of strange attraciors. & ~ r e " 17 shows the Poinear6 are similar to bifurcation d i m a m s in the information thev 0.30. sections of the stranee attractors at F = 0.60 and - - c~=~ provide about the system. Figure 15 is the Lyapunov expo0.50, and 0.70. The Zfferences are readily apparent. nent spectrum corresponding to the bifurcation diagram Conclusion shown in Figure 9, except that to avoid clutter, only the exponents eenerated from IC2 are ~lotted.The zero exwThe damped and driven oscillator in an asymmetric doun e k is notshown. The first occurreke of chaotic behaAor ble-well potential exhibits behavior typical of nonlinear is characterized bv a oositive Lva~unovexwnent at F = dissipative systems as characterized by several methods of 0.365. Regions of the iyapunov exponent spectrum where analysis. The system displays positive Lyapunov expoone exponent is positive correspond to chaotic windows on nents, which indicate chaotic trajectories, for some sets of the bifurcation diagram. Figure 16 is a more detailed parameters. The sum of thes,e exponents measures the Lyapunov exponent spectrum of the region between F = contraction of the phase-space'volume, as calculated from

'

~

~

~

~

~

~

~

..


373

~

~

~

-

Driving Force (F)

~.~~

Fioure exoonents are oldted as a func" - - 15. - The - non-zero - - - Lvaounov -,~- ~ tlon of tne driv ng force ampl'n.de Jslng slanoard parameters and IC2 in the equal ons. The system folows a perloo-00.0 ing roLte lo cnaos, w n c h occbrs at F -0 355, wnere tne largest exponent b e comes greater than zero ~

~

~~

~

~~~

Figure 17. Poincare sections of the strange attractors at F = 0.60, calculated using standard parameters, except (a, green dots) c = 0.30, (b, red dots) c = 0.50, and (c, blue dots) c = 0.70. The fractal dimension of the attractor becomes larger as the damping decreases. This is seen here as an increase in the spacing of the layers as c decreases. ior. The conclusions reached in the first paragraph of this section are valid for all the values of the asymmetry pa! = 0, the Duffing oscillator. It rameter studied, including 3 may, however, change the number of basins of attraction present in phase space. Four were found for P = 0,while only three were located for P = 0.10. This study reflects but a small sampling of the possibilities available for extending the numerical study of nonlinear dynamical systems to the undergraduate level. For this system alone, other techniques, which may have revealed other aspects of the system's behavior, could be used. Variations in other parameters, especially g, may prove interesting. Alarge number of other chaotic systems, suitable for study, have been defined and particularly insightful students may modify or develop techniques or systems.

h i i n g Force (F) Figure 16. Lyapunov exponents plotted versus driving force amplitude in the region F = 0.30 to 0.39. Two sets of exponents. which deoend on the initial conditions. exist in this reaion. Lvaounov exoone'nts generated from IC1 or ld3 are represenied byiquares, while those generated from IC2 are represented by circles. Filled-in symbols denote the most positive exponent. the system's Lie derivative. The period-doubling mute to chaos is followed in accordance with the ~ e i ~ e n b a unumm ber. Strange attractors with fractal dimensions, which increase with decreases in the damping constant, are present in chaotic regions. Asymmetry affects the shape of trajectories; the position of am~litudeiumos. bifurcations. and chaotic and ~eriodic windiws; a d the ihysical form bf strange attractbrs; but it does not affect the overall nature of the system's behav374


Acknowledgment SLT gratefully acknowledges the Sigma Xi Scientific Research Society for support of this project through a Grantin-Aid of Research. CBL would like to acknowledge support from the Idaho State Board of Education (SBOE 92-105). Literature Cited 1. Cleick, J. Chms, Mobingo New S e k n e ; Viking-Penguin:New York, 1987. 2. Rssband, S. T h C h o o t i c D y ~ m i cofNonlimarSysUms; s W8ley:NewYark. 1989. 3. Berge, P.:Pomesu, Y;Vidsl, C. Onbr Within Chaos; Wiley:New York, 1987. 4. Mmn. F Chaotic Vibmtions: Wiley: NewYork. 1987. 5. Baker, C. L.; Gollub, J.P. Chaotic Dynomks, an Inlmductlon: Cambcdge U n i v e ~ sity: NewYork, 1990. 6. msk, H. Diffemntrcll and Difemnrm Eguations thmugh Compu*r Eqerimenta; Springer-Verlsg: NewYork, 1989. 7. De Soura.Maehsdo, s.; Rollins. R. w: Jambe, D. T;nartmao, J. L A m J Phy5. 1990,58(4).321-329. 8. Crebogi, C:Ott, E.;Yorke, J.A. Science 1987.238.632-638. 9. Holmes, P J.,Philos. k n s R Soe London A.1979,292,41948. 10. Scott, S. K,C h m i m l Chaos, hlernorlond SeriesofMonogmphson Ckemislry. 24; Oxford University Press: NewYork, 1991. 11. Ref 5,pp 147-149 and 153-154. 12. Ref. 5,p 149. 13. Wo1f.A.:SuiR, J. B.: Swinney, H. L.;Vaatano,J. A. Physim, 1985,16D,285317. 14. Kaplan, J. L.:Yorke, J. A. In Funetlonol Diffemnllol Equations and UIeAppmximcv fion ofFiredPoints;Peitgen andWalther, Eds.; LectureNofes in Mathematics No. 730; Springer-Verlsg: Berlin, 1978,pp 2W227.

Study of the Chaotic Behavior of a Damped and Driven Oscillator in an

Recommend Documents