Group 1 2012: Difference between revisions

From nldlab
Jump to navigation Jump to search
(Initial import of Latex file via pandoc)
 
m (Protected "Group 1 2012" ([edit=sysop] (indefinite) [move=sysop] (indefinite)))
 
(44 intermediate revisions by one other user not shown)
Line 1: Line 1:
= Objective =
''Group members: Ross Granowski, Aemen Lodhi, Andrew Champion, Suchithra Ravi''


The duffing oscillator is a simple model of a buckling beam which is subjected to lateral periodic forcing. Despite its simplicity, this system exhibits very rich dynamics and is now a classic example of chaos in a necessarily non-autonomous, two-dimensional system. We will attempt to capture the behavior of the system numerically and experimentally via Poincaré sections of trajectories in phase space. We will perform the same analysis with the cosine forcing term replaced by non-sinusoidal periodic functions.
* [[Media:  duffingpreproposal.pdf | Pre-proposal presentation]]
* [[Media:  duffingfinal.pdf | Final presentation]]
* Reports: [[Media: duffingChampion.pdf | Andrew]], [[Media: duffingGranowski.pdf | Ross]]


= Background =
[[File: Tw_duffing.png | frame | Phase portrait of a duffing oscillator [http://en.wikipedia.org/wiki/File:Tw_duffing.png]]]


The Duffing oscillator is modeled by the following nonlinear, second-order, non-autonomous differential equation <math>\ddot{x}+\delta \dot{x}-\alpha x+\beta x^{3}=f \cos(\omega t)</math> for constants <math>\alpha,\beta,\delta,\omega,f>0</math>. Our goal in this experiment is to analyze bifurcations which occur for fixed <math>\alpha,\beta,\delta,\omega</math> as ''f'' is increased. It has been demonstrated numerically that period-doubling bifurcations and chaotic behavior occur in this case as well . However, to the best of our knowledge, this has not been studied experimentally, and nothing has been written about the effect of alternate forcing functions upon the strange attractors in phase space which are observed in the sinusoidal case (although the case of small-amplitude forcing with general periodic functions was studied in ).
<blockquote>We compare a mechanically forced double-well Duffing oscillator with an asymmetric potential to numeric simulations of a symmetric potential system. Analysis is performed with time series data, phase portraits, recurrent Poincare sections, attractor reconstruction from time delay embedding, and Lyapunov exponent estimation from this reconstruction. Our investigation finds that the asymmetric mechanical system can exhibit chaos while for the same forcing parameters symmetric simulation produces limit cycles above or below the chaotic threshold depending on the mechanical potential well used to estimate simulation parameters.
</blockquote>
= Introduction =
 
[[Image:beam_illustration.png|thumb|FIG 0: A common physical interpretation of the double well Duffing oscillator for <math>\alpha > 0</math>|alt=image]]
 
The forced Duffing oscillator is a seminal system for the study of chaotic dynamics and development of analytical and experimental techniques for nonlinear systems. A common physical interpretation for the double-well case of the Duffing oscillator model considered here is a ferromagnetic beam positioned between two magnets while undergoing lateral forcing, as illustrated in FIG 0. This mechanical system elucidates several features of the model Duffing oscillator: a potential for the lateral position of the beam created by the two magnets; non-autonomous forcing resulting in motion of the beam, potentially chaotically between the two wells; and mechanical dissipation within the beam and from air resistance against the beam's motion.
 
Though the Duffing oscillator and the double-well beam system are simple and not of any obvious utility in their own right, they provide a structured model with few, well-defined parameters and clear physical interpretations while still producing chaotic dynamics. This allows intuition into the system's behavior and creates an accessible and familiar setting for experimental design, data analysis, and diagnosis of unexpected results. All these features are valuable both for pedagogical purposes and for the development of new methods and tools for analysis, experimentation, computation, and diagnosis that may be extensible to more complex systems or other scientific or engineering applications.
 
The periodically-forced Duffing oscillator is a system governed by the following nonlinear and nonautonomous second-order ODE (``The Duffing Equation&quot;): <math>\ddot{x}+\delta \dot{x}-\alpha x+\beta x^{3}=f \cos(\omega t)</math> for constants <math>\alpha,\beta,\delta,\omega,f>0</math>. Although it can be viewed as a planar (2 dimensional) system, the nonautonomous periodic term <math>f \cos(\omega t)</math> means that the equation actually determines a smooth vector field on <math>\mathbb{R}^{2}\times \mathbb{S}^{1}</math> (a 3 dimensional space) and hence truly chaotic behavior is possible.


To begin to understand this equation, it is simplest to first analyze the case where <math>\delta=f=0</math>. In this case our equation reduces to <math>\ddot{x}-\alpha x+\beta x^{3}=0</math> which is a Newtonian system with potential <math>V(x)= \frac{\beta }{4}x^{4}-\frac{\alpha}{2} x^{2}</math> and the trajectories are curves of constant energy <math>E(x,\dot{x})=\frac{1}{2}\dot{x}^{2}+V(x)</math>. In this case it is easy to see that we have a saddle at the origin and centers at <math>\left(\pm\sqrt{\frac{\alpha}{\beta}},0\right)</math>. Other than the two homoclinic orbits from the origin to itself and the two centers, we have that all other trajectories are in fact closed orbits.
To begin to understand this equation, it is simplest to first analyze the case where <math>\delta=f=0</math>. In this case our equation reduces to <math>\ddot{x}-\alpha x+\beta x^{3}=0</math> which is a Newtonian system with potential <math>V(x)= \frac{\beta }{4}x^{4}-\frac{\alpha}{2} x^{2}</math> and the trajectories are curves of constant energy <math>E(x,\dot{x})=\frac{1}{2}\dot{x}^{2}+V(x)</math>. In this case it is easy to see that we have a saddle at the origin and centers at <math>\left(\pm\sqrt{\frac{\alpha}{\beta}},0\right)</math>. Other than the two homoclinic orbits from the origin to itself and the two centers, we have that all other trajectories are in fact closed orbits.


Now when we include the dissipation term (i.e. we let <math>\delta>0</math>), the system is again easy to understand: heuristically, trajectories spiral into what are now stable fixed points at <math>\left(\pm\sqrt{\frac{\alpha}{\beta}},0\right)</math>. We can see that in this case <math>\frac{d}{dt}E = \dot{x}(\ddot{x}+\beta x^{3}-\alpha x)=-\delta\dot{x}^{2}\leq 0</math> so the energy decreases (strictly, when <math>\dot{x}\neq 0</math>) along trajectories of the system. The origin is still a saddle in this case. It can be shown that, except for the stable manifolds of the saddle, all trajectories of this system tend toward one of the stable fixed points .
Now when we include the dissipation term (i.e. we let <math>\delta>0</math>), the system is again easy to understand: heuristically, trajectories spiral into what are now stable fixed points at <math>\left(\pm\sqrt{\frac{\alpha}{\beta}},0\right)</math>. We can see that in this case <math>\frac{d}{dt}E = \dot{x}(\ddot{x}+\beta x^{3}-\alpha x)=-\delta\dot{x}^{2}\leq 0</math> so the energy decreases (strictly, when <math>\dot{x}\neq 0</math>) along trajectories of the system. The origin is still a saddle in this case. It can be shown that, except for the stable manifolds of the saddle, all trajectories of this system tend toward one of the stable fixed points.
 
When we let <math>f>0</math>, the behavior of the system becomes quite complicated. For sufficiently small values of ''f'' we find that we have stable limit cycles near what were stable fixed points for <math>f=0</math>, and almost all trajectories eventually converge to one of the two limit cycles (namely all except for the measure zero stable manifolds of the saddle point at the origin). For very large ''f'', it can be shown there is one globally stable limit cycle attracting all trajectories. For intermediate values of ''f'', however, we have chaotic behavior of trajectories in phase space. This is essentially due to the fact that as we increase ''f'' from 0, we eventually cause the stable and unstable manifolds of the saddle at <math>(0,0)</math> to intersect, and hence to intersect infinitely many times. The presence of infinitely many homoclinic intersection points in (a compact region of) phase space causes the trajectories to have very complex topological structures. These can be seen most clearly via the Poincaré section (sampling the orbit at time intervals of length <math>\frac{2\pi}{\omega}</math>), which displays intricatre fractal ``strange attractors&quot; for these chaotic values of ''f''.
 
In our experiment, we constructed a mechanical system whose behavior is approximately governed by the above equation. In essence, we used two magnets placed symmetrically about the tip of a strip of shim stock to give us a system behaving like <math>\ddot{x}=-\delta \dot{x}+\alpha x-\beta x^{3}</math>. Here ''x'' denotes the displacement of the tip of the stock, with the vertical beam giving us ''x''=0 and the magnets being located at <math>\pm\sqrt{\frac{\alpha}{\beta}}</math>. We can think of the tip of the stock as a particle trapped in a ``double well&quot; potential. We then used a motor to generate the periodic forcing by shaking the apparatus in the lateral direction.
 
Our apparatus was modeled on that of Moon and Holmes <ref name="MoonHolmes"/>, where they use essentially the same setup (although inverted - which we are assuming makes almost no difference, since gravity never appears in any of the ``governing equations&quot;). In the Moon and Holmes paper, they show that the Duffing equation is indeed a good model for this system via a clever reduction (using Galerkin&rsquo;s method) of a system of PDEs. This analysis, although interesting, is somewhat complicated and does not significantly improve our intuition for the system, so we refer the interested reader to the Moon and Holmes paper for further details. Indeed, the readers seeking a more detailed mathematical analysis of this system should conslult the papers of Holmes <ref name="Holmes"/> and Moon and Holmes<ref name="MoonHolmes"/>.
 
Our goal at the outset was to analyze the bifurcations of the physical system in the <math>(f,\omega)</math> plane for fixed values <math>\alpha,\beta,\delta</math>. We state here that the Duffing equation is already nondimensionalized. Later, we will use our experimental results to approximate the values of <math>\alpha, \beta, \delta</math> appearing in the eqation. We the run computer simulations with these parameters and compare them to our data.
 
= Methods =


The orbits become much more complicated for <math>f>0</math>. For small ''f'', the behavior is well-understood (see  for details). In this case, the system closely resembles the non-forced system: the stable fixed points become attracting orbits of period <math>\frac{2\pi}{\omega}</math> and we still have something like a saddle near the origin. For large ''f'', the behavior is also well-understood, since in this case the forcing term ``dominates&quot; the other terms. In between these two extremes, the system becomes extremely complicated. In order to understand the behavior of the system with a ``mid-range&quot; non-zero forcing term, <math>f>0</math>, its easier to study and analyze the Poincaré section of the trajectories.
== Experimental Setup ==


== Poincaré section of the oscillator ==
[[Image:nldExperimentSchematic.png|thumb|FIG 1: Components and Connection of Experimental Apparatus|alt=image]]


The Poincaré section is obtained by plotting <math>(x(t),y(t))</math> whenever ''t'' is a multiple of <math>\frac{2\pi}{\omega}</math>. This is the same thing as sampling the system at the same phase in every cycle. The Poincaré section for this system shows how the points fall on a fractal set and is very sensitive to initial conditions indicating the possibility of sustained chaos in the experiment. The Poincaré sections at low ''f'' is topologically equivalent to the sections at <math>f=0.</math> P0 has 3 fixed points and the manifolds do not intersect. As ''f'' increases, however, the manifolds become tangential and can intersect transversely. The creation of such a global homoclinic bifurcation leads to the presence of infinite transversal homoclinic orbits inidicating the presence of infinite number of unstable periodic orbits (of arbitrary periods). It also indicates the presence of bounded nonperiodic motion. All these together create the structure of the strange attractor. Of course, this also leads to the theoretical question of whether the strange attractor is simply an artifact of noise or does it really exist. The question of existence of strange attractors is a very pertinent one, but is not being explored here and is only mentioned for the sake of completeness.
For the mechanical system itself (the source of the data), we attached a slender metallic beam to a rigid frame made from 80/20. The beam was cut from shim stock of thickness .007 inches (.0001778 m). Although we had originally planned to use electromagnets, we would have had to weight the end of the beam in order to produce the desired buckling behavior. Weighting the stock correctly proved rather difficult, so insead we used very powerful stationary magnets and an unweighted beam. Two rare earth magnets were attached symmetrically to the top of the frame. The strength of the magnets was such that the metallic beam always buckled in the direction of one of magnets, with either configuration (being buckled to the left or to the right) constituting a stable equilibrium. The position of the beam in the center is an unstable equilibrium. Since the focus of this experiment was on understanding the forced vibrations of the beam, we then applied an external force to the entire framework with a motor whose position was controlled through the computer.


== Analog with the double well ==
We note here that for controlling the motion of the motor we specified the frequency (in Hz &ndash; generally 3, 4, or 5 Hz &ndash; although we wrote a program in LabView which allowed us to generate motion at arbitrary frequencies), the forcing function (a sine wave), and the amplitude of oscillation. The amplitude of oscillation was measured in terms of ``counts,&quot; where 1000 counts = .006 m.


This system is analogous to a damped particle in a double well when the entire system is shaken. The two wells represent the two stable states of the buckled beam and the hump represents the saddle point when the beam is straight. When we apply a small force to one of the stable points, minor oscillations around the bottom of the well are produced. If a large force is applied, the particle is simply thrown from one well to the other. In the intermediate case (i.e., for forcing values in between) the particle may start out at the hump and chaotically move between the wells, because the particle falls one way or the other based on the timing and value of the forcing function as well as the particle&rsquo;s energy.
To measure the displacement of the tip, we decided to use strain gauges, which were attached near base of beam. The strain gauges &ndash; essentially a type of variable resistor &ndash; were used to create a voltage differential which was then amplified and sent to the computer where it was plotted as a function of time. Our sampling rate here was 10,000 Hz (i.e. we recorded 10,000 voltage measurements per second). The strain gauges are basically measuring the local deflection at the base of the beam (in the immediate vicinity of the strain gauge). Thus we worried that the resulting voltage measurements would not be adequate, since we are really just concerned with the position of the tip of the beam. We had originally inteded to use a webcam as well, but we found that the strain gauges actually had good accuracy and precision, minimal drift, and adequate range. Hence a camera was not needed. Measuring the tip position versus voltage reading, we found an approximate linear correspondence between voltage and distance: 1.4V corresponds to approximately .045 m.


== Research questions ==
A conceptual schematic of the experimental setup is presented in FIG. 1, above.


The scientific question we are trying to answer is whether different forcing functions produce qualitatively different chaotic behavior and bifurcations. Does the geometry of the strange attractor change based on the type of forcing function we use? Does the chaotic behavior persist for all kinds of forcing functions? These questions are important because existing literature has primarily focused on single cosine functions, but since this system is nonlinear, even a simple sum of cosines might lead to qualitatively different results. Real world signals typically contain noise or may have other forms (e.g., square or triangular). So, generalization of studies is possible only after we understand whether trajectories are completely different when the forcing functions are changed.
== Evolution of the Apparatus ==


= Experimental methods =
This apparatus went through several iterations as we attempted to correct various problems and improve the design.


We use the following experimental setup to explore the dynamics of the Duffing oscillator.
For the first iteration (see FIG. 2) we had a wheeled cart on a low-friction track (borrowed from the undergraduate Physics I lab). The cart carried the 80/20 frame and we held the motor stationary and attached the motor arm to the cart. There measurement setup was fine, but there were several mechanical problems with this setup: the motor was difficult to attach to the cart and there was a problematic amount of friction coming from the wheels of the cart (the weight of the frame on top was causing a wheel to ``stick&quot;).


== Electro-mechanical system concept ==
For the second iteration (see FIG 3) we replaced the cart and track with an air track. This time we held the motor arm stationary while attaching the frame to the top of the motor. This was a significant improvement over the first iteration (basicallly eliminating the friction problem). But we still had several problems: most importantly, the shim stock would strike the frame when its motion became chaotic. For large amplitudes and/or frequencies of oscillation, this setup would also cause intense vibrations of the table and the entire apparatus.


We attach a slender metallic beam to a rigid framework. Two electromagnets are attached at the base of the framework. The strength of the magnets is such that the metallic beam is pulled in the direction of one of magnets. Either configuration is stable. The position of the beam in the center is an unstable equilibrium as small disturbances to the framework cause the beam to buckle to one of the magnets. Since the goal of this experiment is to understand the forced vibration of the beam we apply an external force to the entire framework. Figure {fig:Experimental-setup-concept} gives the concept of the experimental setup.
For the third and final iteration (see FIG 4) we addressed these problems by widening the 80/20 frame to accomodate chaotic motion and we mounted the apparatus on cinderblocks which were adhered to the floor and to one another. We still have problems with this setup, however. There are problems with effectively cooling the motor and the CAN amplifier. The motor&rsquo;s current limitiations are also preventing us from driving the system hard enough to observe sustained chaotic motion. There is also the issue of asymmetry in the system due to a difference in magnet strength and the natural bend in the stock, neither of which have we been able to effectively eliminate or compensate for. This is apparent in our data below.


<blockquote><blockquote>[[Image:experimentalSetup|frame|none|alt=image]]
<div class="center" style="width:auto; margin-left:auto; margin-right:auto;">
</blockquote>
{|
{Experimental setup concept{fig:Experimental-setup-concept}}
| [[Image:nldFirstIteration.jpg|thumb|FIG 2: First Iteration of the Setup|alt=image]]
</blockquote>
| [[Image:nldSecondIteration.jpg|thumb|FIG 3: Second Iteration of the Setup|alt=image]]
== Inputs and outputs ==
| [[Image:nldThirdIteration.png|thumb|FIG 4: Third Iteration of the Setup|alt=image]]
|}
</div>
 
= Results =
 
===Early Apparatus Video===
<videoflash>9Lt8kfEg5FA</videoflash>
===Beam Becoming Stuck in a Well===
<videoflash>7gS2Z92mufg</videoflash>
===Final Apparatus Video===
<videoflash>NCeBxAi-D6E</videoflash>
 
== Parameter Estimation ==
 
Following Berger and Nunes<ref name="BergerNunes"/>, we estimate the values of <math>\alpha, \beta,\text{ and }\delta</math> appearing in the Duffing equation using our experimental data. We will then use these estimated values (together with the known values of ''f'' and <math>\omega</math>) in a computer simulation, which we compare to the results of a particular trial.
 
According to our data, we have fixed points at approximately <math>\pm\sqrt{\frac{\alpha}{\beta}}\approx \pm .7 V</math> (the strain gauges would shift over time, but this is the position once we normalize by subtracting the mean from the data). Writing <math>x(t)=-\sqrt{\frac{\alpha}{\beta}}+\phi(t)</math>, substituting this into the Duffing equation and ignoring the forcing term, we find that to first order in <math>\phi</math> we have <math>\ddot{\phi}(t)=-\delta\dot{\phi}(t)-2\alpha\phi(t)</math>. We know that this gives us damped simple harmonic motion with frequency <math>\omega_{1}=\sqrt{2\alpha(1-\delta^{2}/8\alpha)}</math>. Solving in terms of <math>\delta</math> gives us <math>\delta =\sqrt{8\alpha(1-\omega_{1}^{2}/2\alpha)}</math>. From our experimental data we get an approximate value for <math>\omega_{1}</math> (based on the period of oscillation in the small amplitude case) and guess <math>\alpha</math> such that <math>\omega_{0}=\sqrt{2\alpha}</math> is only slightly larger than <math>\omega_{1}</math>.
 
From a set of data where the stock remained oscillating with low amplitude in the right well (time series pictured in FIG. 5 above), we found that <math>\omega_{1,\text{left}}=2.0\text{cycles/second}</math>. Thus we have <math>\beta=\frac{\alpha}{.49}</math> and <math>\delta= \sqrt{8\alpha(1-2/\alpha)}</math>. We observed chaotic behavior for forcing with amplitude 2100 count (so approximately .392 V, using .006 m per count and 1.4V per .045m) and frequency 4Hz so <math>f=\text{amplitude}*\text{frequency}^{2}\approx 6.272</math>. Setting <math>\alpha=2.01</math> is physically reasonable, since then <math>\omega_{0}=\sqrt{2\alpha}\approx 2.004</math> is the undamped frequency (this should be just a little bigger than the damped frequency <math>\omega_{1,\text{left}}</math>). For the left well we find <math>\omega_{1_\text{right}}=4.76\text{cycles/sec}</math>. Again <math>\beta=\frac{\alpha}{.49}</math> and now <math>\delta= \sqrt{8\alpha(1-22.6576/2\alpha)}</math>. In this case (for our simulation) we set <math>\alpha=11.33</math> (which gives us <math>\omega_{0}\approx 4.7602</math>).Averaging the two values of <math>\omega_{1}</math> we get <math>\omega_{1, \text{avg}}=3.38</math>. Again <math>\beta=\frac{\alpha}{.49}</math> and now <math>\delta= \sqrt{8\alpha(1-11.424/2\alpha)}</math>. Now for our simulation with these parameters we set <math>\alpha=5.714</math> (which gives us <math>\omega_{0}\approx 3.3805</math>). {Plots}
 
In what follows we present several plots from one of our runs in which we observed chaotic behavior. The plot in FIG. 6 is a time series plot of the voltage measurements collected during the run. One magnetic well is located at around <math>-.4 V</math> (on the ''y''-axis) and the other is at around <math>-1.8 V</math>. One can see that the motions of the tip are highly irregular and indeed appear random. We then passed this data through a Butterworth filter to clean up the high-frequency noise in the signal. The resulting data was then numerically differentiated (we simply used the difference quotients <math>\frac{f(x+h)-f(x)}{h}</math> since our time steps <math>h</math> were extremely small, namely ten-thousandths of a second). Then we plotted the smoothed signal versus its derivative to obtain the picture of the trajectory in phase space (<math>x\text{ vs. }\dot{x}</math>), which is shown in FIG. 7. The Poincaré section in FIG. 8 is obtained from the data plotted in FIG. 7 by plotting one point per oscillation of the motor. The motor was oscillating at 4 Hz and we were sampling at 10,000 Hz so we plotted only the first of every 2,500 data points. The points were then colored to indicate the sequence in which they were plotted.
 
<div class="center" style="width:auto; margin-left:auto; margin-right:auto;">
{|
| [[Image:1500count_3Hz_leftWell_LimitCycle.png|400px|thumb|FIG 5: A ``small f&quot; limit cycle, ''y''-axis units are Volts, ''x''-axis units are (<math>\frac{1}{10,000}</math>)ths of a second|alt=image]]
| [[Image:2100count_4Hz_TimeSeries.png|400px|thumb|FIG 6: Chaotic behavior for sinusoidal forcing of amplitude .0126 m at 4Hz. ''y''-axis units are Volts, ''x''-axis units are (<math>\frac{1}{10,000}</math>)ths of a second|alt=image]]
|-
| [[Image:2100count_4Hz_PhasePlot.png|400px|thumb|FIG 7: Chaotic behavior for sinusoidal forcing of amplitude .0126 m at 4Hz. ''y''-axis units are Volts per second, ''x''-axis units are Volts|alt=image]]
| [[Image:2100count_4Hz_PoincareColored.png|400px|thumb|FIG 8: Poincaré section of trajectory in FIG. 7. Sinusoidal forcing of amplitude .0126 m at 4Hz. ''y''-axis units are Volts per second, ''x''-axis units are Volts. The color indicates the order in which the points were plotted.|alt=image]]
|}
</div>
 
== Attractor Reconstruction. Estimating the Largest Lyapunov Exponent. ==
 
We followed the method outlined in Fraser and Swinney<ref name="FraserSwinney"/> for reconstructing attractors from time series. We use the time series <math>X(t)</math> from our chaotic trial whose graph is shown in FIG. 6 and we plot <math>t\mapsto (X(t), X(t+\tau), X(t+2\tau))</math>, where <math>\tau</math> is chosen to be the first local minimum of the mutual information (see FIG. 9) for <math>(X(t), X(t+\tau))</math>. We do not have the space here to explain this concept in detail, but, heuristically, we want <math>X(t)</math> and <math>X(t+\tau)</math> to be as statistically independent as possible. The 3-dimensional embedding of the reconstructed attractor which results is shown in FIG. 10.
 
Using the reconstructed attractor plot, we then followed the method in Rosenstein et. al.<ref name="Rosenstein"/> for estimating the largest Lyapunov exponent. We use an algorithm to find the nearest neighbor of each point on the reconstructed attractor, breaking our time series into a set of closest pairs of points. Then letting <math>d_{j}(i)</math> denote the distance between the <math>j</math>th pair of nearest neighbors at time <math>i\Delta t</math>, we plot the line <math>y(i)=\frac{1}{\Delta t}\langle \ln d_{j}\rangle</math>, where <math>\langle\cdot \rangle</math> is the average over ''j''. The largest Lyapunov exponent is then given by the slope of the least squares fit to this line.
 
With this method, we found the largest Lyapunov exponent to be .0196 - which is positive, as we would expect for a system exhibiting chaotic behavior (see FIG. 11).
 
<div class="center" style="width:auto; margin-left:auto; margin-right:auto;">
{|
| [[Image:miMinimization.png|400px|thumb|FIG 9: Minimizing the Mutual Information|alt=image]]
| [[Image:3dReconstructedAttractor.png|400px|thumb|FIG 10: The reconstructed attractor. The units for each axis are Volts.|alt=image]]
|-
| [[Image:lyapunovExponent.png|400px|thumb|Plot of <math>\frac{1}{\Delta t}\langle \ln d_{j}(i) \rangle</math> (on the ''y''-axis) versus <math>i</math> (on the ''x''-axis). The slope of the least-squares fit line gives us an estimate of the largest Lyapunov exponent.|alt=image]]
| [[Image:2dReconstructedAttractor.png|400px|thumb|FIG 10b: The 2D-embedded reconstructed attractor. The units for each axis are Volts.|alt=image]]
|}
</div>


The following parameters will be provided as inputs to the experimental setup:
= Discussion =


# Periodic forcing function frequency <math>\omega</math>
The phase plots of the simulations run with the parameters estimated in III.B are shown in FIG. 12, FIG. 13, and FIG. 14 below.
# Magnitude of the force <math>F</math>
# Strength of the electromagnet <math>E</math>


Superpositions of forcing functions and non-sinusoidal forcing functions, such as saw waves, will also be explored. The following signals will be measured as output from the experimental setup:
Compare these plots to that of the one we were attempting to simulate (shown in FIG. 7) and you will notice several significant discrepancies. For one thing, only the simulation using the ``right well&quot; parameters showed any significant chaotic behavior, and even in that case the chaotic transient ended fairly quickly in a limit cycle. None of these plots had enough chaotic behavior to generate any meaningful Poincaré section or to use the method of Lyapunov exponent estimation via attractor reconstruction.


# The displacement of the metallic beam from the center position as a function of time <math>x(t)</math>
The other fact that needs to be pointed out is that the experimental plot in FIG. 7 is extremely asymmetric (the trajectories spending a disproportionately long time in the right well, which is apparently much larger than the left well). The simulation plots, on the other hand, are quite symmetric. We believe this asymmetry is due to the natural bend of the shim stock in the direction of the magnet corresponding to the right-hand potential well. There are also noticeable differences in the strength of the magnets, which certainly contributes to this observed asymmetry as well.


== Physical setup ==
<div class="center" style="width:auto; margin-left:auto; margin-right:auto;">
{|
| [[Image:simulation_leftWell_params_phasePlot.png|400px|thumb|FIG 12: The plot in phase space of the simulation starting from <math>(-.7,0)</math>, using the parameters obtained from the ``small ''f''&quot; oscillations in the right well.|alt=image]]
| [[Image:simulation_rightWell_params_phasePlot.png|400px|thumb|FIG 13: The plot in phase space of the simulation starting from <math>(-.7,0)</math>, using the parameters obtained from the ``small ''f''&quot; oscillations in the left well.|alt=image]]
|-
| [[Image:simulation_averageWell_params_phasePlot.png|400px|thumb|FIG 14: The plot in phase space of the simulation starting from <math>(-.7,0)</math>, using the parameters obtained from the average of ``small ''f''&quot; oscillations in the left and right wells.|alt=image]]
|
|}
</div>


The physical setup of the experiment will consist of the following equipment, whose logical arrangement is shown in Figure {fig:Experiment-physical-setup}:
= Conclusion =


# Air track and cart to act as the rigid framework shown in Figure {fig:Experimental-setup-concept}
We must acknowledge that the forced Duffing oscillator is now something of a classical example of a chaotic dynamical system, and it has been thoroughly studied in theory and experiment. Although we did not necessarily do anything new, we found it to be a very deep and complicated system. It was a great chance to apply the theory which was presented in the class. In the course of our experiment and analysis we learned about a number of interesting and fruitful techniques for the study of chaotic dynamics.
# Metallic beam, composed of shim stock with brass weights attached at the end as necessary to observe the desired dynamics
# Electromagnets
# Linear motor
# Signal generator for the forcing function
# Power supply
# Strain gauge to measure the displacement of the beam from the center. We will employ a strain gauge as well as a camera to record displacement in preliminary experiments to determine tradeoff in accuracy, sensor drift and computational complexity between the two methods. Visual measurements may require a large amount of computational storage space for long running experiments. Using a strain gauge may allow measurement of the curvature of the metallic beam while requiring less data storage and analysis. The strain gauge will be attached to a computer to record the displacement of the beam&rsquo;s end from the center.
# A plexiglass setup to house the entire experimental setup (optional). We will measure the effect of air resistance on the experiment and if the effect is non-negligible we will house the entire experimental setup in a plexiglass case to isolate experiment.


<blockquote><blockquote>[[Image:expsetup.pdf|frame|none|alt=image]]
We still have number of questions regarding the Duffing oscillator and our experiment which we find very interesting. We would like to explore several different directions in the future. In particular, we would like to analyze forced oscillators with asymmetric potential wells<ref name="Asymmetric"/>, which can be modeled by setting <math>V(x)= \frac{\beta }{4}x^{4}+\frac{\epsilon_{1}}{3}x^{3}-\frac{\alpha}{2} x^{2}+\epsilon_{2}x</math>. We would also like to study the case of nonsinusoidal periodic forcing<ref name="VariousPeriodicForces"/>, where we replace <math>\cos</math> with some other function <math>\rho(\omega t)</math> which is possibly piecewise smooth (or merely piecewise <math>C^{\alpha}(\mathbb{R})</math>). Exploring either of these questions will require further refining the experimental setup, continued study of numerical simulations, and some new mathematical/theoretical techniques.
</blockquote>
{Physical components and electrical connections of the experimental apparatus} {fig:Experiment-physical-setup}
</blockquote>
Other factors such as air resistance, background magnetic fields, and loss of calibration of the strain gauge may impact experimental results and require modification of the experiment or sensor hardware.


== Cost of experiment ==
= References =


We do not expect any costly procurements for this experiment. We expect that the equipment required to carry out this experiment will be available in any basic physics laboratory.
<references>
<ref name="BergerNunes">Berger, J. and G. Nunes Jr. (1997). A mechanical duffing oscillator for the undergraduate laboratory. ''Am. J. Phys.''. 65, 275-296.</ref>
<ref name="FraserSwinney">Fraser, A. and H. Swinney (1986). Independent coordinates for strange attractors from mutual information. ''Phys. Rev. A,'' 33, 1134-1140.</ref>
<ref name="Holmes">Holmes, P. (1979). A nonlinear oscillator with a strange attractor. ''Phil. Trans. R. Soc. Lond''. A 292, 419–448.</ref>
<ref name="MoonHolmes">Moon, F. and P. Holmes (1979). A magnetoelastic strange attractor. ''Journal of Sound and Vibration''. 65, 841-846.</ref>
<ref name="VariousPeriodicForces">Ravichandran, V., V. Chinnathambi, and S. Rajasekar (2006). Effect of various periodic forces on duffing oscillator. ''Pramana 67''(2), 351–356.</ref>
<ref name="Rosenstein">Rosenstein, M.T. and Collins, J.J. and De Luca, C.J. (1993). A practical method for calculating largest Lyapunov exponents from small data sets. "Physica D: Nonlinear Phenomena". 65, 117-134.</ref>
<ref name="Asymmetric">S. Jeyakumari, V. Chinnathambi, S. Rajasekar, M.A.F. Sanjuan (2011). Vibrational resonance in an asymmetric duffing oscillator. "I. J. Bifurcation and Chaos". 21, 275-286.</ref>
</references>

Latest revision as of 22:50, 28 October 2014

Group members: Ross Granowski, Aemen Lodhi, Andrew Champion, Suchithra Ravi

Phase portrait of a duffing oscillator [1]

We compare a mechanically forced double-well Duffing oscillator with an asymmetric potential to numeric simulations of a symmetric potential system. Analysis is performed with time series data, phase portraits, recurrent Poincare sections, attractor reconstruction from time delay embedding, and Lyapunov exponent estimation from this reconstruction. Our investigation finds that the asymmetric mechanical system can exhibit chaos while for the same forcing parameters symmetric simulation produces limit cycles above or below the chaotic threshold depending on the mechanical potential well used to estimate simulation parameters.

Introduction

image
FIG 0: A common physical interpretation of the double well Duffing oscillator for <math>\alpha > 0</math>

The forced Duffing oscillator is a seminal system for the study of chaotic dynamics and development of analytical and experimental techniques for nonlinear systems. A common physical interpretation for the double-well case of the Duffing oscillator model considered here is a ferromagnetic beam positioned between two magnets while undergoing lateral forcing, as illustrated in FIG 0. This mechanical system elucidates several features of the model Duffing oscillator: a potential for the lateral position of the beam created by the two magnets; non-autonomous forcing resulting in motion of the beam, potentially chaotically between the two wells; and mechanical dissipation within the beam and from air resistance against the beam's motion.

Though the Duffing oscillator and the double-well beam system are simple and not of any obvious utility in their own right, they provide a structured model with few, well-defined parameters and clear physical interpretations while still producing chaotic dynamics. This allows intuition into the system's behavior and creates an accessible and familiar setting for experimental design, data analysis, and diagnosis of unexpected results. All these features are valuable both for pedagogical purposes and for the development of new methods and tools for analysis, experimentation, computation, and diagnosis that may be extensible to more complex systems or other scientific or engineering applications.

The periodically-forced Duffing oscillator is a system governed by the following nonlinear and nonautonomous second-order ODE (``The Duffing Equation"): <math>\ddot{x}+\delta \dot{x}-\alpha x+\beta x^{3}=f \cos(\omega t)</math> for constants <math>\alpha,\beta,\delta,\omega,f>0</math>. Although it can be viewed as a planar (2 dimensional) system, the nonautonomous periodic term <math>f \cos(\omega t)</math> means that the equation actually determines a smooth vector field on <math>\mathbb{R}^{2}\times \mathbb{S}^{1}</math> (a 3 dimensional space) and hence truly chaotic behavior is possible.

To begin to understand this equation, it is simplest to first analyze the case where <math>\delta=f=0</math>. In this case our equation reduces to <math>\ddot{x}-\alpha x+\beta x^{3}=0</math> which is a Newtonian system with potential <math>V(x)= \frac{\beta }{4}x^{4}-\frac{\alpha}{2} x^{2}</math> and the trajectories are curves of constant energy <math>E(x,\dot{x})=\frac{1}{2}\dot{x}^{2}+V(x)</math>. In this case it is easy to see that we have a saddle at the origin and centers at <math>\left(\pm\sqrt{\frac{\alpha}{\beta}},0\right)</math>. Other than the two homoclinic orbits from the origin to itself and the two centers, we have that all other trajectories are in fact closed orbits.

Now when we include the dissipation term (i.e. we let <math>\delta>0</math>), the system is again easy to understand: heuristically, trajectories spiral into what are now stable fixed points at <math>\left(\pm\sqrt{\frac{\alpha}{\beta}},0\right)</math>. We can see that in this case <math>\frac{d}{dt}E = \dot{x}(\ddot{x}+\beta x^{3}-\alpha x)=-\delta\dot{x}^{2}\leq 0</math> so the energy decreases (strictly, when <math>\dot{x}\neq 0</math>) along trajectories of the system. The origin is still a saddle in this case. It can be shown that, except for the stable manifolds of the saddle, all trajectories of this system tend toward one of the stable fixed points.

When we let <math>f>0</math>, the behavior of the system becomes quite complicated. For sufficiently small values of f we find that we have stable limit cycles near what were stable fixed points for <math>f=0</math>, and almost all trajectories eventually converge to one of the two limit cycles (namely all except for the measure zero stable manifolds of the saddle point at the origin). For very large f, it can be shown there is one globally stable limit cycle attracting all trajectories. For intermediate values of f, however, we have chaotic behavior of trajectories in phase space. This is essentially due to the fact that as we increase f from 0, we eventually cause the stable and unstable manifolds of the saddle at <math>(0,0)</math> to intersect, and hence to intersect infinitely many times. The presence of infinitely many homoclinic intersection points in (a compact region of) phase space causes the trajectories to have very complex topological structures. These can be seen most clearly via the Poincaré section (sampling the orbit at time intervals of length <math>\frac{2\pi}{\omega}</math>), which displays intricatre fractal ``strange attractors" for these chaotic values of f.

In our experiment, we constructed a mechanical system whose behavior is approximately governed by the above equation. In essence, we used two magnets placed symmetrically about the tip of a strip of shim stock to give us a system behaving like <math>\ddot{x}=-\delta \dot{x}+\alpha x-\beta x^{3}</math>. Here x denotes the displacement of the tip of the stock, with the vertical beam giving us x=0 and the magnets being located at <math>\pm\sqrt{\frac{\alpha}{\beta}}</math>. We can think of the tip of the stock as a particle trapped in a ``double well" potential. We then used a motor to generate the periodic forcing by shaking the apparatus in the lateral direction.

Our apparatus was modeled on that of Moon and Holmes <ref name="MoonHolmes"/>, where they use essentially the same setup (although inverted - which we are assuming makes almost no difference, since gravity never appears in any of the ``governing equations"). In the Moon and Holmes paper, they show that the Duffing equation is indeed a good model for this system via a clever reduction (using Galerkin’s method) of a system of PDEs. This analysis, although interesting, is somewhat complicated and does not significantly improve our intuition for the system, so we refer the interested reader to the Moon and Holmes paper for further details. Indeed, the readers seeking a more detailed mathematical analysis of this system should conslult the papers of Holmes <ref name="Holmes"/> and Moon and Holmes<ref name="MoonHolmes"/>.

Our goal at the outset was to analyze the bifurcations of the physical system in the <math>(f,\omega)</math> plane for fixed values <math>\alpha,\beta,\delta</math>. We state here that the Duffing equation is already nondimensionalized. Later, we will use our experimental results to approximate the values of <math>\alpha, \beta, \delta</math> appearing in the eqation. We the run computer simulations with these parameters and compare them to our data.

Methods

Experimental Setup

image
FIG 1: Components and Connection of Experimental Apparatus

For the mechanical system itself (the source of the data), we attached a slender metallic beam to a rigid frame made from 80/20. The beam was cut from shim stock of thickness .007 inches (.0001778 m). Although we had originally planned to use electromagnets, we would have had to weight the end of the beam in order to produce the desired buckling behavior. Weighting the stock correctly proved rather difficult, so insead we used very powerful stationary magnets and an unweighted beam. Two rare earth magnets were attached symmetrically to the top of the frame. The strength of the magnets was such that the metallic beam always buckled in the direction of one of magnets, with either configuration (being buckled to the left or to the right) constituting a stable equilibrium. The position of the beam in the center is an unstable equilibrium. Since the focus of this experiment was on understanding the forced vibrations of the beam, we then applied an external force to the entire framework with a motor whose position was controlled through the computer.

We note here that for controlling the motion of the motor we specified the frequency (in Hz – generally 3, 4, or 5 Hz – although we wrote a program in LabView which allowed us to generate motion at arbitrary frequencies), the forcing function (a sine wave), and the amplitude of oscillation. The amplitude of oscillation was measured in terms of ``counts," where 1000 counts = .006 m.

To measure the displacement of the tip, we decided to use strain gauges, which were attached near base of beam. The strain gauges – essentially a type of variable resistor – were used to create a voltage differential which was then amplified and sent to the computer where it was plotted as a function of time. Our sampling rate here was 10,000 Hz (i.e. we recorded 10,000 voltage measurements per second). The strain gauges are basically measuring the local deflection at the base of the beam (in the immediate vicinity of the strain gauge). Thus we worried that the resulting voltage measurements would not be adequate, since we are really just concerned with the position of the tip of the beam. We had originally inteded to use a webcam as well, but we found that the strain gauges actually had good accuracy and precision, minimal drift, and adequate range. Hence a camera was not needed. Measuring the tip position versus voltage reading, we found an approximate linear correspondence between voltage and distance: 1.4V corresponds to approximately .045 m.

A conceptual schematic of the experimental setup is presented in FIG. 1, above.

Evolution of the Apparatus

This apparatus went through several iterations as we attempted to correct various problems and improve the design.

For the first iteration (see FIG. 2) we had a wheeled cart on a low-friction track (borrowed from the undergraduate Physics I lab). The cart carried the 80/20 frame and we held the motor stationary and attached the motor arm to the cart. There measurement setup was fine, but there were several mechanical problems with this setup: the motor was difficult to attach to the cart and there was a problematic amount of friction coming from the wheels of the cart (the weight of the frame on top was causing a wheel to ``stick").

For the second iteration (see FIG 3) we replaced the cart and track with an air track. This time we held the motor arm stationary while attaching the frame to the top of the motor. This was a significant improvement over the first iteration (basicallly eliminating the friction problem). But we still had several problems: most importantly, the shim stock would strike the frame when its motion became chaotic. For large amplitudes and/or frequencies of oscillation, this setup would also cause intense vibrations of the table and the entire apparatus.

For the third and final iteration (see FIG 4) we addressed these problems by widening the 80/20 frame to accomodate chaotic motion and we mounted the apparatus on cinderblocks which were adhered to the floor and to one another. We still have problems with this setup, however. There are problems with effectively cooling the motor and the CAN amplifier. The motor’s current limitiations are also preventing us from driving the system hard enough to observe sustained chaotic motion. There is also the issue of asymmetry in the system due to a difference in magnet strength and the natural bend in the stock, neither of which have we been able to effectively eliminate or compensate for. This is apparent in our data below.

image
FIG 2: First Iteration of the Setup
image
FIG 3: Second Iteration of the Setup
image
FIG 4: Third Iteration of the Setup

Results

Early Apparatus Video

<videoflash>9Lt8kfEg5FA</videoflash>

Beam Becoming Stuck in a Well

<videoflash>7gS2Z92mufg</videoflash>

Final Apparatus Video

<videoflash>NCeBxAi-D6E</videoflash>

Parameter Estimation

Following Berger and Nunes<ref name="BergerNunes"/>, we estimate the values of <math>\alpha, \beta,\text{ and }\delta</math> appearing in the Duffing equation using our experimental data. We will then use these estimated values (together with the known values of f and <math>\omega</math>) in a computer simulation, which we compare to the results of a particular trial.

According to our data, we have fixed points at approximately <math>\pm\sqrt{\frac{\alpha}{\beta}}\approx \pm .7 V</math> (the strain gauges would shift over time, but this is the position once we normalize by subtracting the mean from the data). Writing <math>x(t)=-\sqrt{\frac{\alpha}{\beta}}+\phi(t)</math>, substituting this into the Duffing equation and ignoring the forcing term, we find that to first order in <math>\phi</math> we have <math>\ddot{\phi}(t)=-\delta\dot{\phi}(t)-2\alpha\phi(t)</math>. We know that this gives us damped simple harmonic motion with frequency <math>\omega_{1}=\sqrt{2\alpha(1-\delta^{2}/8\alpha)}</math>. Solving in terms of <math>\delta</math> gives us <math>\delta =\sqrt{8\alpha(1-\omega_{1}^{2}/2\alpha)}</math>. From our experimental data we get an approximate value for <math>\omega_{1}</math> (based on the period of oscillation in the small amplitude case) and guess <math>\alpha</math> such that <math>\omega_{0}=\sqrt{2\alpha}</math> is only slightly larger than <math>\omega_{1}</math>.

From a set of data where the stock remained oscillating with low amplitude in the right well (time series pictured in FIG. 5 above), we found that <math>\omega_{1,\text{left}}=2.0\text{cycles/second}</math>. Thus we have <math>\beta=\frac{\alpha}{.49}</math> and <math>\delta= \sqrt{8\alpha(1-2/\alpha)}</math>. We observed chaotic behavior for forcing with amplitude 2100 count (so approximately .392 V, using .006 m per count and 1.4V per .045m) and frequency 4Hz so <math>f=\text{amplitude}*\text{frequency}^{2}\approx 6.272</math>. Setting <math>\alpha=2.01</math> is physically reasonable, since then <math>\omega_{0}=\sqrt{2\alpha}\approx 2.004</math> is the undamped frequency (this should be just a little bigger than the damped frequency <math>\omega_{1,\text{left}}</math>). For the left well we find <math>\omega_{1_\text{right}}=4.76\text{cycles/sec}</math>. Again <math>\beta=\frac{\alpha}{.49}</math> and now <math>\delta= \sqrt{8\alpha(1-22.6576/2\alpha)}</math>. In this case (for our simulation) we set <math>\alpha=11.33</math> (which gives us <math>\omega_{0}\approx 4.7602</math>).Averaging the two values of <math>\omega_{1}</math> we get <math>\omega_{1, \text{avg}}=3.38</math>. Again <math>\beta=\frac{\alpha}{.49}</math> and now <math>\delta= \sqrt{8\alpha(1-11.424/2\alpha)}</math>. Now for our simulation with these parameters we set <math>\alpha=5.714</math> (which gives us <math>\omega_{0}\approx 3.3805</math>). {Plots}

In what follows we present several plots from one of our runs in which we observed chaotic behavior. The plot in FIG. 6 is a time series plot of the voltage measurements collected during the run. One magnetic well is located at around <math>-.4 V</math> (on the y-axis) and the other is at around <math>-1.8 V</math>. One can see that the motions of the tip are highly irregular and indeed appear random. We then passed this data through a Butterworth filter to clean up the high-frequency noise in the signal. The resulting data was then numerically differentiated (we simply used the difference quotients <math>\frac{f(x+h)-f(x)}{h}</math> since our time steps <math>h</math> were extremely small, namely ten-thousandths of a second). Then we plotted the smoothed signal versus its derivative to obtain the picture of the trajectory in phase space (<math>x\text{ vs. }\dot{x}</math>), which is shown in FIG. 7. The Poincaré section in FIG. 8 is obtained from the data plotted in FIG. 7 by plotting one point per oscillation of the motor. The motor was oscillating at 4 Hz and we were sampling at 10,000 Hz so we plotted only the first of every 2,500 data points. The points were then colored to indicate the sequence in which they were plotted.

image
FIG 5: A ``small f" limit cycle, y-axis units are Volts, x-axis units are (<math>\frac{1}{10,000}</math>)ths of a second
image
FIG 6: Chaotic behavior for sinusoidal forcing of amplitude .0126 m at 4Hz. y-axis units are Volts, x-axis units are (<math>\frac{1}{10,000}</math>)ths of a second
image
FIG 7: Chaotic behavior for sinusoidal forcing of amplitude .0126 m at 4Hz. y-axis units are Volts per second, x-axis units are Volts
image
FIG 8: Poincaré section of trajectory in FIG. 7. Sinusoidal forcing of amplitude .0126 m at 4Hz. y-axis units are Volts per second, x-axis units are Volts. The color indicates the order in which the points were plotted.

Attractor Reconstruction. Estimating the Largest Lyapunov Exponent.

We followed the method outlined in Fraser and Swinney<ref name="FraserSwinney"/> for reconstructing attractors from time series. We use the time series <math>X(t)</math> from our chaotic trial whose graph is shown in FIG. 6 and we plot <math>t\mapsto (X(t), X(t+\tau), X(t+2\tau))</math>, where <math>\tau</math> is chosen to be the first local minimum of the mutual information (see FIG. 9) for <math>(X(t), X(t+\tau))</math>. We do not have the space here to explain this concept in detail, but, heuristically, we want <math>X(t)</math> and <math>X(t+\tau)</math> to be as statistically independent as possible. The 3-dimensional embedding of the reconstructed attractor which results is shown in FIG. 10.

Using the reconstructed attractor plot, we then followed the method in Rosenstein et. al.<ref name="Rosenstein"/> for estimating the largest Lyapunov exponent. We use an algorithm to find the nearest neighbor of each point on the reconstructed attractor, breaking our time series into a set of closest pairs of points. Then letting <math>d_{j}(i)</math> denote the distance between the <math>j</math>th pair of nearest neighbors at time <math>i\Delta t</math>, we plot the line <math>y(i)=\frac{1}{\Delta t}\langle \ln d_{j}\rangle</math>, where <math>\langle\cdot \rangle</math> is the average over j. The largest Lyapunov exponent is then given by the slope of the least squares fit to this line.

With this method, we found the largest Lyapunov exponent to be .0196 - which is positive, as we would expect for a system exhibiting chaotic behavior (see FIG. 11).

image
FIG 9: Minimizing the Mutual Information
image
FIG 10: The reconstructed attractor. The units for each axis are Volts.
image
Plot of <math>\frac{1}{\Delta t}\langle \ln d_{j}(i) \rangle</math> (on the y-axis) versus <math>i</math> (on the x-axis). The slope of the least-squares fit line gives us an estimate of the largest Lyapunov exponent.
image
FIG 10b: The 2D-embedded reconstructed attractor. The units for each axis are Volts.

Discussion

The phase plots of the simulations run with the parameters estimated in III.B are shown in FIG. 12, FIG. 13, and FIG. 14 below.

Compare these plots to that of the one we were attempting to simulate (shown in FIG. 7) and you will notice several significant discrepancies. For one thing, only the simulation using the ``right well" parameters showed any significant chaotic behavior, and even in that case the chaotic transient ended fairly quickly in a limit cycle. None of these plots had enough chaotic behavior to generate any meaningful Poincaré section or to use the method of Lyapunov exponent estimation via attractor reconstruction.

The other fact that needs to be pointed out is that the experimental plot in FIG. 7 is extremely asymmetric (the trajectories spending a disproportionately long time in the right well, which is apparently much larger than the left well). The simulation plots, on the other hand, are quite symmetric. We believe this asymmetry is due to the natural bend of the shim stock in the direction of the magnet corresponding to the right-hand potential well. There are also noticeable differences in the strength of the magnets, which certainly contributes to this observed asymmetry as well.

image
FIG 12: The plot in phase space of the simulation starting from <math>(-.7,0)</math>, using the parameters obtained from the ``small f" oscillations in the right well.
image
FIG 13: The plot in phase space of the simulation starting from <math>(-.7,0)</math>, using the parameters obtained from the ``small f" oscillations in the left well.
image
FIG 14: The plot in phase space of the simulation starting from <math>(-.7,0)</math>, using the parameters obtained from the average of ``small f" oscillations in the left and right wells.

Conclusion

We must acknowledge that the forced Duffing oscillator is now something of a classical example of a chaotic dynamical system, and it has been thoroughly studied in theory and experiment. Although we did not necessarily do anything new, we found it to be a very deep and complicated system. It was a great chance to apply the theory which was presented in the class. In the course of our experiment and analysis we learned about a number of interesting and fruitful techniques for the study of chaotic dynamics.

We still have number of questions regarding the Duffing oscillator and our experiment which we find very interesting. We would like to explore several different directions in the future. In particular, we would like to analyze forced oscillators with asymmetric potential wells<ref name="Asymmetric"/>, which can be modeled by setting <math>V(x)= \frac{\beta }{4}x^{4}+\frac{\epsilon_{1}}{3}x^{3}-\frac{\alpha}{2} x^{2}+\epsilon_{2}x</math>. We would also like to study the case of nonsinusoidal periodic forcing<ref name="VariousPeriodicForces"/>, where we replace <math>\cos</math> with some other function <math>\rho(\omega t)</math> which is possibly piecewise smooth (or merely piecewise <math>C^{\alpha}(\mathbb{R})</math>). Exploring either of these questions will require further refining the experimental setup, continued study of numerical simulations, and some new mathematical/theoretical techniques.

References

<references> <ref name="BergerNunes">Berger, J. and G. Nunes Jr. (1997). A mechanical duffing oscillator for the undergraduate laboratory. Am. J. Phys.. 65, 275-296.</ref> <ref name="FraserSwinney">Fraser, A. and H. Swinney (1986). Independent coordinates for strange attractors from mutual information. Phys. Rev. A, 33, 1134-1140.</ref> <ref name="Holmes">Holmes, P. (1979). A nonlinear oscillator with a strange attractor. Phil. Trans. R. Soc. Lond. A 292, 419–448.</ref> <ref name="MoonHolmes">Moon, F. and P. Holmes (1979). A magnetoelastic strange attractor. Journal of Sound and Vibration. 65, 841-846.</ref> <ref name="VariousPeriodicForces">Ravichandran, V., V. Chinnathambi, and S. Rajasekar (2006). Effect of various periodic forces on duffing oscillator. Pramana 67(2), 351–356.</ref> <ref name="Rosenstein">Rosenstein, M.T. and Collins, J.J. and De Luca, C.J. (1993). A practical method for calculating largest Lyapunov exponents from small data sets. "Physica D: Nonlinear Phenomena". 65, 117-134.</ref> <ref name="Asymmetric">S. Jeyakumari, V. Chinnathambi, S. Rajasekar, M.A.F. Sanjuan (2011). Vibrational resonance in an asymmetric duffing oscillator. "I. J. Bifurcation and Chaos". 21, 275-286.</ref> </references>