
  
\documentclass[11pt]{article}
\usepackage{html}
\pagestyle{empty}
  
% Fit to page specifications:
%
%       Defaults for article:   textwidth = 5 in
%                               textheight = 7.25 in
%                               hoffset = 0 in
%                               voffset = 0 in
%
%       These settings center the text on the page, with right and
%       left margins of 1.75 in and top and bottom margins of 2 in
%       for standard U.S. letter size paper.
%       Changing hoffset moves the text to the *right* by hoffset.
%       Changing voffset moves the text *down* by voffset.
  
\textwidth 6.5in
\textheight 9.5in
  
\begin{document}
  
\begin{rawhtml}
        <center><img src=../../gif_files/GRAD_BLUE_LINE_f.gif></center>
\end{rawhtml}
 
\begin{center}
        {\LARGE\bf The Duffing Oscillator }
\end{center}
 
\begin{rawhtml}
        <center><img src=../../gif_files/GRAD_BLUE_LINE_f.gif></center>
\end{rawhtml}


The original Duffing oscillator (Duffing, 1918) was introduced 
in relation to the single (spatial) modes of vibration of 
a steel beam subjected 
to external periodic forces. It has now become one of
the standard prototype of forced systems thanks to
its simplicity and yet the richness of its solution.
This is specially true of the modified model with
a double well potential.

\begin{rawhtml}
<center>
<img src=../temp.gif width="460" height="262">
</center>
\end{rawhtml}


The model is
described in terms of a potential well
$$
V(x) = \beta x^2 / 2 + \alpha x^4 / 4 
$$
$\beta$ defines the shape of the potential. It is usually
chosen to be +1 or -1. $\alpha$ is 
a measure of the non-linearity. The latter is usually
taken to be 1. The steel beam is assumed to oscillate
harmonically for small displacement $x$, nonlinearity setting
in at larger $x$ values. For $\beta = 1$, the potential 
exhibits one minimum at $x = 0.0$.  For $\beta = -1$, this
 potential is a symmetric two well potential corresponding
to the physics of a steel beam influenced by two magnets
as in the figure above.


From this potential, the following force is
derived: 
$F_{beam}(x) = - \frac{ dV } { dx } = 
             - \beta x - \alpha x^3  $.
The steel beam is also subjected to a dissipation
term, $F_{dissipation} = - \mu v$, proportional to
the velocity. A driving term is also applied,
$F_{driving} =  A cos( \omega t )$, which corresponds
to a back-and-forth motion of the point of suspension
of the beam.

The model is summarized in the two following ODEs:
$$
\frac{dx}{dt} = v
$$ 
$$
\frac{dv}{dt} = \left( - \beta x - \alpha x^3 - \mu v + A cos( \omega t ) \right) / m 
$$
The parameters of the model are $\alpha, \beta, \mu, \omega,$ and $A$.
The mass is $m$. 
The time is represented by $t$.


References abound for the Duffing Oscillator. For instance,
you can look online at 
\begin{rawhtml}
<a href=http://www.mcasco.com/pattr1.html>Double Well Duffing</a>
\end{rawhtml}. 
Google shows thousands of links.

You can also read {\em Nonlinear Dynamics: A Two-way Trip from Physics to Math}
by Hernan G. Solari, Mario A. Natiello, and Gabriel B. Mindlin, 
published by The Institute of Physics,  September, 1996.
Dr. Solari and Mindlin were at Drexel in the past as a 
researcher and graduate student respectively.

An understanding of the Duffing Oscillator
 in term of {\em topology of the orbits} can be
found in {\em The topology of Chaos}, Robert Gilmore and
Marc Lefrand, Wiley, 2002.

\subsection{Coding a Solution of the Duffing Oscillator}

The Duffing Oscillator system is described in terms of 
 a set of two simple coupled
ODEs to solve.  This solution can be obtained suffciently
accurately via the Runge Kutta 4th order
method.

The program to implement the numerical solution of
the Duffing Oscillator
can be derived from the Runge Kutta 4th order program
to solve the mass on a spring problem. 

Adapt this code according to the following steps

\begin{itemize}

\item Modify the RHS of the differential equation and the total energy
         according to the {\em Duffing Oscillator} equations

\item Set initial conditions to
         $x_0 = 1.0$ and $v_0 = 0.0$ in {\em set\_initial\_conditions()}

\item Define the model parameters as 
      $\beta = 1$, $\alpha, \mu, A$ as zero, and $\omega = 1 $
      in {\em define\_model()}

\item Use $dt = 0.02$ and $t_{max} = 25.0$ in {\em define\_time\_grid()}

\item Code a piece of logic to implement the recognition of three line
    arguments: -p, -i, and -t n

    They have the following functions:
   \begin{description}
     \item[-p] Print trajectory data. Default not to print.
     \item[-i] Read in the initial conditions $x_0$ and $v_0$ from
              {\em <stdin> or cin} until $EOF$ is reached.
     \item[-t n] To set $t_{max} = n 2 \pi / \omega $, $n$ times
                the period of the external driving  term 
   \end{description}

\item Implement the functionality of the three line arguments above.

    The -i option is implemented via an infinite $for( ~~; ~~; ~~ )$ loop
    with a $ if( ~~~~ ) break; $ statement 

\end{itemize}


\subsection{Harmonic Oscillator Limit}

{\em \bf Harmonic Oscillator limit -- $\beta = 1$ case}

With $\beta = 1$ and $\alpha, \mu, A$ all equal to zero 
 this system
is a simple harmonic oscillator (mass=1).

The solution of the model produces a phase space 
portrait in $v(t)$ versus $x(t)$ of simple
concentric circles surrounding the {\em Fixed Point} at 
the origin. The {\em Fixed Point} of a set of ODEs 
 is at coordinates that render the
RHS of the ODEs equal to zero.

\begin{rawhtml}
<img src=../hand_pictures/HO_phase_space.gif height=200 width=200 border=0>
\end{rawhtml}


\begin{description}
\item[Exercise \#1:]
This system is 
a harmonic oscillator. Draw the phase space portrait of this
solution. You should write a {\em HO setup code} to feed initial
conditions in the ODE solver code. The program
should sprinkle initial conditions along the positive $x$ axis.
The latter must be written with 
a {\em program control} ( -i ) to control whether or not to read
in initial conditions from <stdin>.

Use the general modular RK4 ODE solver program
to solve the Duffing system, using $\alpha =0.0$, $\beta = 1$,
$\mu = 0.0$, $A = 0.0$ and $\omega = 1$ in the code.  
Use $dt = 0.02$. 

Check that the energy of each trajectory is conserved by the
numerical solution.
\end{description}


\bigskip

{\em \bf Dissipative System}

With $\mu$
not zero, the system becomes a {\em damped harmonic
oscillator}. This term takes energy away from the 
motion of the steel beam. All trajectories spiral
into the attractive {\em Fixed Point} at the origin.

\begin{rawhtml}
<img src=../hand_pictures/HO_dissipative_v_x.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#2:]
Draw the phase space portrait of this
solution. Use the {\em HO setup code} to feed initial
conditions in the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha = 0.0$, $\beta = 1$,
$\mu = 0.1$, $A = 0.0$ and $\omega = 1$ in the code. 

Plot the energy of the system along these trajectories
as functions of time.

Note the energy dissipation.

\end{description}



\subsection{Modified Harmonic Oscillator - One Well System}

Setting $\alpha = 1.0$ (usual value in standard analysis of
the model) changes the potential well shape. 

\begin{rawhtml}
<img src=../hand_pictures/HO_NL.gif height=200 width=300 border=0>
\end{rawhtml}

This changes 
the shape of the phase space portrait. The extra
stiffness of the beam {\em squares}  the trajectories by imparting
a stronger impact on the beam to push it in 
opposite direction than the force resulting from
a parabolic potential by itself could have.

\begin{rawhtml}
<img src=../hand_pictures/NL_phase_space.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#3:]
Reproduce the 
phase space portrait as obtained for the HO in
the previous exercises.
Use the {\em HO setup code} to specify the
initial conditions. 
Use the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha = 1.0$, $\beta = 1$,
$\mu = 0.0$, $A = 0.0$ and $\omega = 1$ in the code.
\end{description}


\bigskip

{\em \bf Dissipative System}

The phase space portrait with $\mu$ not zero
 shows that the origin, 
$x = 0$, $v = 0$ remains an attractive 
{\em equilibrium point}
or {\em fixed point}. This point
makes the RHS of the ODEs equal to zero. It corresponds
to a minimum energy point. Trajectories spiral toward this 
{\em fixed point}.

\begin{rawhtml}
<img src=../hand_pictures/NL_dissipative_v_x.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#4:]
 Reproduce the phase space portrait above by
using initial conditions generated by the {\em HO setup code}.
Use $\alpha = 1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 0.0$ and $\omega = 1$ in the 
code implementing the RK4 solution of {\em Duffing Oscillator}.

Plot the energy corresponding to each of the initial
conditions used in building the phase space portrait
above. Note the energy dissipation.
\end{description}

 


\subsection{Two Wells System}

Setting $\beta=-1$ dramatically 
changes the shape of the potential.
It becomes a {\em Two Wells Potential}.

\begin{rawhtml}
<img src=../hand_pictures/2w_potential.gif height=180 width=300 border=0>
\end{rawhtml}

This potential harbors two minima at $x$ = +1 and -1 ( $\alpha = 1.0$ ) and
a maximum at $x = 0$. The phase space portrait of the solution 
of this problem is very different than that of the single well
problem.

\begin{rawhtml}
<img src=../hand_pictures/2w_v_x.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#5:]
Write a {\em 2 wells setup code} to generate
initial conditions to launch trajectories that will illustrate
the behavior of the solutions of the {\em Two Wells Potential}.
These initial conditions will be fed into the ODE solving code.
These initial conditions should be sprinkled along
the $x$ axis between the wells and to the right of the wells.

Use $\alpha = 1.0$, $\beta = -1$,
$\mu = 0.0$, $A = 0.0$ and $\omega = 1$ in the 
code implementing the RK4 solution of {\em Duffing Oscillator}.

Check that the energy is conserved, i.e., constant, for each trajectory.
\end{description}



\subsection{Homoclinic Orbit}

The {\em Fixed Point} at the origin is unstable 
(maximum of potential). It
becomes interesting to study the trajectories
initiating near the origin. 

There is a special trajectory going through
the origin called {\em homoclinic orbit}
or {\em separatrix}. It separates the 
inner orbits circling the stable
{\em Fixed Points} (minimum of potential)
from the trajectories circling both
wells. 

\begin{rawhtml}
<img src=../hand_pictures/2w_homoclinic.gif height=200 width=200 border=0>
\end{rawhtml}

The directions in which the {\em homoclinic orbit} spreads
out from the origin can be found by analyzing 
trajectories based on initial
conditions distributed on a small circle centered on the origin.
The trajectories based on initial conditions near the 
unstable manifolds will be dragged outward while those
based on initial conditions near the 
stable manifolds will come in toward the unstable point
before being dragged outward.

\begin{rawhtml}
<img src=../hand_pictures/distortion_homoclinic.gif height=200 width=200 border=0>
\end{rawhtml}

By lengthening the integration time the different trajectories
approach  the {\em homoclinic orbit}.


\begin{description}
\item[Exercise \#6:]
Write a code to generate 100 initial conditions on
a small circle around the origin to feed into the
ODE solver code. Integrate the ODEs to a small 
$t_{max}$ until larger values to illustrate
the deformation of the circle of initial
conditions.

Use the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha = 1.0$, $\beta = -1$,
$\mu = 0.0$, $A = 0.0$ and $\omega = 1$ in the code.

\item[Exercise \#7:]
Run the integration over larger yet $t_{max}$ to see
the shape of the {\em homoclinic orbit}.

\end{description}





\bigskip

{\em \bf Dissipative System}

Dissipation drives the trajectories in spirals 
toward either of the two stable {\em Fixed Points}
 located at the minimum of the potential.
The homoclinic orbit also spirals toward both
stable {\em Fixed Points}.

\begin{rawhtml}
<img src=../hand_pictures/2w_dissipation_v_x.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#8:]
Use the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha = 1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 0.0$ and $\omega = 1$ in the code.

Feed the initial conditions generated via the
{\em 2 wells setup code} in the ODE solver code.

Plot a phase space portrait of the solutions.


\item[Exercise \#9:]
Superpose the homoclinic orbit generated via
        initial conditions on a circle onto the graph
         above
\end{description} 




\subsection{Driven One Well Duffing Oscillator}

Driving the Duffing Oscillator with either $\beta = 1$ or
$\beta = -1$  produces incredibly complex behavior, leading
to chaos for some range of the parameters.

At low driving strength the one well Duffing Oscillator
 exhibits attractive 
{\em Limit Cycles}. This is a trajectory in phase space
toward which trajectories converge following a transient
period. The region over which initial conditions lead to 
trajectories converging to the {\em Limit Cycle} is the 
{\em Basin of Attraction}. The following is a sketch 
of the attractive {\em Limit Cycle} for $A = 2.0$.

\begin{rawhtml}
<img src=../hand_pictures/1w_driven_2.0.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#10:]
Draw the phase space portrait of this
solution. Use the {\em HO setup code} to feed initial
conditions in the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha =1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 2.0$ and $\omega = 1$ in the code. 

Note that the different trajectories seem to converge 
toward a common trajectory, a {\em Limit Cycle}


\item[Exercise \#11:]
Refine the graph of the {\em Limit Cycle}. To do so, integrate
the ODEs starting from an arbitrary initial conditions,
i.e., $(x_0,v_0) = ( 1.0,0.0),$ for a large $t_{max}$
via the option -t 800, i.e., 800 periods of the external field. \\

Restart the integration from suitable initial conditions
for another long integration time interval to refine the {\em Limit Cycle}
calculation.

Do it again, i.e., run the ODE integrator again!


Plot $x(t)$ and $v(t)$ as functions of $t$.   Note that the period
of the {\em Limit Cycle} is that of the external forcing term.

\end{description}


\bigskip

{\em \bf Larger Driving Strengths}

Larger driving strengths can produce very different
solutions in the driven one well Duffing Oscillator. For instance
the following sketch illustrates the two 
{\em Limit Cycles}  at a driving
strength $A=4.0$. 

\begin{rawhtml}
<img src=../hand_pictures/1w_driven_4.0.gif height=200 width=200 border=0>
\end{rawhtml}

These two {\em Limit Cycles} are mirror images of each other.
If a solution of a lesser symmetry than that of the potential
exists, a solution of opposite symmetry should also exist.

\begin{description}
\item[Exercise \#12:] 
Draw the phase space portrait of this
solution. Use the {\em HO setup code} to feed initial
conditions in the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha =1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 4.0$ and $\omega = 1$ in the code. 

Note that the different trajectories seem to converge 
toward two possible trajectories, two {\em Limit Cycles}

\item[Exercise \#13:]
Define better these {\em Limit Cycles}. To do so, integrate
the ODEs starting from an arbitrary initial conditions,
i.e., $(x_0,v_0) = ( 1.0,0.0),$ for a large $t_{max}$
via the option -t 800.

Restart the integration from suitable initial conditions
for another long integration time interval to refine the 
two {\em Limit Cycles} calculations.

\item[Exercise \#14:]
  Repeat the calculations above for other initial
  conditions in the range $(x_0,v_0) = ( 0.9,0.0)$
  down to $(0.3,0.0) $.
  
  It is clear that each {\em Limit Cycle} has its own 
  {\em Basin of Attraction}.
\end{description}

The following sketch illustrates the {\em Limit Cycle}
generated by the Duffing Oscillator driven with 
a strength $A=6.0$.

\begin{rawhtml}
<img src=../hand_pictures/1w_driven_6.0.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#15:]
Draw the phase space portrait of this
solution.  Use $\alpha =1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 6.0$ and $\omega = 1$ as parameters.
\end{description}




\subsection{Chaotic Solution of the Driven One Well Duffing Oscillator}

Driving harder the One Well Duffing Oscillator eventually
results in the system going chaotic. A characteristic
of this solution is that a chaotic trajectory will
fill the allowed region of phase space. The following is
a sketch of a trajectory in phase space  for
the driven Duffing Oscillator with driving 
strength $A = 35.5$.

\begin{rawhtml}
<img src=../hand_pictures/1w_chaos_v_x_35.5.gif height=200 width=200 border=0>
\end{rawhtml}

It is clear that such solution is very different that we
saw previously!

\begin{description}
\item[Exercise \#16:]
Draw the phase space portrait of this
solution. Use the {\em HO setup code} to feed initial
conditions in the general modular RK4 ODE solver program
to solve the Duffing system. Use $\alpha =1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 35.5$ and $\omega = 1$ in the code.
Increase the integration time ( -t 200 option) and watch 
for the phase space to begin to fill.

 
\item[Exercise \#17:]
Draw the phase space again, this time starting from 
 $(x_0,v_0) = ( 1.0,0.0)$ ( or whatever point in phase space!) 
and integrating
over long integration times ( -t 2000 for instance ). Phase
space fills again!

\end{description}


\bigskip

{\em \bf Order within Chaos}

The solution above is chaotic.  Plotting $x(t)$ and $v(t)$ 
versus $t$ shows that the solution has lost any trace
of periodic behavior (despite the fact that the driving
term is perfectly periodic!). The patterns do not
repeat themselves.

The solution ventures all over the phase space but 
the motion remains within a finite region of space.
The solution exhibits a very sensitive dependence over
initial conditions - this is a hallmark of chaotic behavior.   

\begin{description}
\item[Exercise \#18:]
Solve the Duffing Oscillator ( $\alpha = 1.0$, $\beta = 1$,
$\mu = 0.1$, $A = 35.5$ and $\omega = 1$) starting from 
 $(x_0,v_0) = ( 1.0,0.0)$  and integrating
over long integration times ( -t 2000 for instance ).
Restart the solution starting from 
$( x_0, v_0 ) = ( 1.00001, 0.0 )$. Plot $ x(t) $ versus
$t$ together for each case over increasing time intervals
and watch how the two solutions diverge from
each other (exponential divergence).

\end{description}

{\em \bf Stroboscopic Map - Strange Attractor }

All this points to the chaotic (random) aspects of the 
solution. To see {\em order} within this solution 
requires another way of thinking. We need to look at 
the solution in a phase space that includes the time $t$ axis,
i.e., $x(t), v(t), t$.

\begin{description}
\item[Exercise \#19:]
Plot the chaotic solution in $x(t), v(t), t$ space
( gnuplot subcommand splot ).
\end{description}

The system is driven by an external driving term of
angular frequency $\omega = 2 \pi f = 1.0$, corresponding
to a period $T = 2 \pi / \omega$.  Imagine asking 
where is the steel beam  and how fast does it go 
at each period of the external driving term
  and recording this info in a phase space like graph
$v(n T)$ versus $x(n T)$.  Each time $t=nT$ defines
a plane in $v(n T)$ versus $x(n T)$ which the trajectory
traverses at {\em one point}. So a trajectory
will generate {\em one dot} per time interval
$T$.  This is called a {\em Stroboscopic Surface of Section} 
or a {\em Stroboscopic Map}.

A {\em Limit cycle} would generate very few
distinguishable points in such a construct since
the trajectory repeats a pattern. The chaotic
solution of the Duffing Oscillator generates
the following picture.

\begin{rawhtml}
<img src=../hand_pictures/1w_driven_stroboscopic.gif height=200 width=200 border=0>
\end{rawhtml}

This is a {\em  Strange Attractor}.

\begin{description}
\item[Exercise \#20:]
Write a {\em Stroboscopic Code} which accepts the
output of the general modular RK4 ODE solver program
to solve the Duffing system and outputs the values
of $x(t)$, $v(t)$ and $energy(t)$ at $t = n T$.
To do so find when the time has crossed  $t = n T$
and perform a linear interpolation to find exact
values of $x(t)$, $v(t)$ and $energy(t)$ at
the crossing.

\item[Exercise \#21:]
Draw the image above with 20,000 points.  
\end{description}










\subsection{Driven Two Well Duffing Oscillator}

The behavior of the Driven Two Well Duffing Oscillator
is even more stunning that the Driven One Well system.


\bigskip

{\em \bf Limit Cycles}

At small driving amplitudes the system stays bounded in either
of the two wells and exhibits {\em Limit Cycles}.
The following sketch shows the solution at driving strength $A = 0.22$.
Note that the amplitude of the driving term
used here is smaller than in the
one well system since a new measure of strength is set by
the barrier height between the two wells.

\begin{rawhtml}
<img src=../hand_pictures/2w_forced_0.22.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#22:]
Reproduce the phase space portrait above by
running the general modular RK4 ODE solver program
to solve the Duffing system over long period of time
( -t 800 ) for the parameters  $\alpha =1.0$, $\beta = -1$,
$\mu = 0.25$, $A = 0.22$  and $\omega = 1$ starting at initial
conditions $(x_0,v_0) = ( 1.0,0.0)$.

\item[Exercise \#23:]
Refine the image of the {\em Limit Cycle}.  
\end{description}





\bigskip

{\em \bf Transient Chaos}

At larger driving amplitudes the system may traverse
from one well to the other. At not such a high
 driving strength, the systems exhibits temporarily
a chaotic behavior before settling in a  {\em Limit Cycle}.
This is called {\em transient chaos}.

The following sketch shows the solution at driving strength $A = 0.245$.

\begin{rawhtml}
<img src=../hand_pictures/2w_forced_0.245.gif height=200 width=200 border=0>
\end{rawhtml}

\begin{description}
\item[Exercise \#24:]
Reproduce the phase space portrait above by
running the general modular RK4 ODE solver program
to solve the Duffing system over long period of time
( -t 800 ) for the parameters  $\alpha =1.0$, $\beta = -1$,
$\mu = 0.25$, $A = 0.245$  and $\omega = 1$.

\item[Exercise \#25:]
Plot $x(t)$ versus $t$ in the transient time period
to illustrate the oscillations from one well to the other,
a characteristic of the chaotic behavior.

\item[Exercise \#26:]
Refine the image of the {\em Limit Cycle}.  
\end{description}




\bigskip

{\em \bf Chaos}

At larger yet driving amplitudes the system has an easier time 
traversing
from one minimum to the other. The system may then become chaotic.
A characteristic of the chaoticity will then be for the steel
beam to bounce from
one well to another in an erratic manner.

As in the case of the one well system, the solution now fills
in the phase space.

\begin{rawhtml}
<img src=../hand_pictures/2w_forced_0.4.gif height=200 width=200 border=0>
\end{rawhtml}

 Only a stroboscopic surface of section can reveal the order in this solution.


\begin{description}
\item[Exercise \#27:]
Reproduce the stroboscopic surface of section  above by
running the general modular RK4 ODE solver program
to solve the Duffing system over long period of time
( -t 2000 ) for the parameters  $\alpha =1.0$, $\beta = -1$,
$\mu = 0.25$, $A = 0.4$  and $\omega = 1$ and piping the results
in the {\em stroboscopic code}.
Use initial conditions $x_0 = 1.0$ and $v_0 = 0.0$.

\item[Exercise \#28:]
Plot $x(t)$ versus $t$ over small  time periods
to illustrate the oscillations from one well to the other,
a characteristic of the chaotic behavior.

\item[Exercise \#29:]
Refine the image of the stroboscopic surface of section. 

\item[Exercise \#30:]
Start from different initial conditions; what happens to
the figure formed in  the stroboscopic surface of section.
\end{description}





\subsection{Strange Attractor}

The stroboscopic surface of section  we produced in
the previous section is a {\em Strange Attractor}.

It is an {\em attractor} in the sense that trajectories
originating from various initial conditions converge
to it. It is {\em Strange} in the sense that it is neither 
a continuous line in $x$ - $v$ phase space (like a 
{\em Limit Cycle} which forms a continuous trajectory
in phase space) nor is it filling the $x$ - $v$ plane.
It is a {\em Fractal}.

As in the case of the Driven One Well chaotic 
solution, this solution exhibits an extreme
sensitivity to initial conditions. This is
the hallmark of chaos.

This sensitivity to initial conditions and the
attractiveness of the {\em Strange Attractor} 
can be illustrated
by following a bunch of trajectories based on
initial conditions that are fairly close together. 
These trajectories should all converge to the
the  {\em Strange Attractor} if observed in
a stroboscopic map at $t = nT$.   Therefore they will
build the shape of the {\em Strange Attractor}
point by point. These points  will be far apart
on the {\em Strange Attractor} due to the 
exponential divergence of the orbits.
This  {\em Strange Attractor} is really {\em Strange}!

\begin{description}
\item[Exercise \#31:]
Write a {\em 2 wells chaos initial setup code} to generate
initial conditions to launch trajectories that scan a 
grid $N_x \times N_v$ on a small rectangle
in $x$-$v$ space. Use $N_x = 50 $, $N_v = 50 $ spanning
$x_0 = [-1.1,-1.0]$ and $v_0 = [-0.1,0.1]$ to start with.
These initial conditions will be fed into the ODE solving code.
Use $\alpha = 1.0$, $\beta = -1$,
$\mu = 0.25$, $A = 0.4$ and $\omega = 1$ in the 
code implementing the RK4 solution of {\em Duffing Oscillator}.
Make sure that the code outputs the last $x$ and $v$ values 
when the trajectory data are not printed out ( by not 
using the -p option ).

Run the code for 1, 2, 3, 4, 5, 10, 20, 100 outside driving
term periods ( via the option -t n ) thereby producing
stroboscopic maps data to be plotted separately.  This will 
illustrates divergence of the trajectories and the attraction of the
{\em Strange Attractor}.

\item[Exercise \#32:]
Repeat the step above with more initial conditions,
e.g., 20000, to refine the image of the {\em Strange Attractor}.

\item[Exercise \#33:]
Start from a different region of phase space.

\end{description}





\subsection{Helper Codes and Model Parameters}



{\bf \em Codes}

\begin{rawhtml}

<ol>
 
<li> The Runge-Kutta 4th order modular ODE solver for the
mass on a spring problem
<br><br>
 
<li> Initial conditions generator code,
<a href=../beta_plus/HO_setup.c>HO_setup.c</a> or
<a href=../beta_plus/HO_setup.cpp>HO_setup.cpp</a> ,
to display Harmonic Oscillator and one-well modified phase space
<br><br>
 
<li> Initial conditions generator code, 
<a href=../beta_minus/2w_setup.c>2w_setup.c </a> or
<a href=../beta_minus/2w_setup.cpp>2w_setup.cpp </a> ,
 to display
phase space of two-well potential phase space
<br><br>
 
<li> Initial conditions generator code, 
<a href=../beta_minus/2w_setup_homoclinic.c>2w_setup_homclinic.c</a> or
<a href=../beta_minus/2w_setup_homoclinic.cpp>2w_setup_homclinic.cpp</a> ,
to display homoclinic orbit in two-well potential
<br><br>
 
<li> Initial conditions generator code,
<a href=../beta_minus/2w_setup_chaotic_attractor.c>2w_setup_chaotic_attractor.c</a> or
<a href=../beta_minus/2w_setup_chaotic_attractor.cpp>2w_setup_chaotic_attractor.cpp</a> ,
to display the chaotic attractor in a dissipative driven two-well potential
<br><br>

<li> Stroboscopic analysis code,
  <a href=../stroboscopic.c>stroboscopic.c</a> or
  <a href=../stroboscopic.cpp>stroboscopic.cpp</a> ,
to display the chaotic attractor in a dissipative driven one- or
   two-well potential. It assumes omega = 1.0
<br><br>

 
</ol>


ODE solver code: <a href=../duffing.c>duffing.c</a> or 
<a href=../duffing.cpp>duffing.cpp</a>

\end{rawhtml}

\bigskip

{\bf \em Model Parameters}

The following are suggested parameters that lead to interesting results.

One-well Duffing Oscillator:

Parameters: $\alpha = 1.0$, $\beta = 1.0$, $\mu = 0.1$, $\omega = 1.0$ 
and $A$ in the range [2.0, 35.5]


Two-well Duffing Oscillator

Parameters: $\alpha = 1.0$, $\beta = -1.0$, $\mu = 0.25$, $\omega = 1.0$ 
and $A$ with values
 [0.22, 0.245, 0.4]


\end{document}


