DOING PHYSICS WITH PYTHON

CHAOTIC DYNAMICAL SYSTEMS

A DRIVEN DAMPED PENDULUM

Ian Cooper

matlabvisualphysics@gmail.com

DOWNLOAD DIRECTORIES FOR PYTHON CODE

       Google drive

       GitHub

cs_006_01.py

Solves the differential equation that governs the dynamics of a simple damped and driven pendulum. The nonlinear ODE for the pendulum system is extremely sensitive to the model’s parameters and initial conditions. All input parameters much be chosen with care. Values mostly used in the simulations are those used by J.R. Taylor in his excellent book Classical Physics.

cs_006_02.py

Lyapunov exponent

cs_006_03.py

Bifurcation diagram

cs_006_04.py

Poincare section plot: for different values of the drive strength , you will need to change the X and Y limits and the ranges for the plotting  and  in some of the plots.

CHAOTIC DYNAMICAL SYSTEMS

Chaotic phenomena appear in many real-life situations from meteorology, medicine, biological systems, ecological systems, fluid dynamics, economics, and many other fields. This article will consider a sinusoidally driven and damped simple pendulum (DDP) as a dynamical system that exhibits chaotic motion.

It is important to realize that chaotic behaviour and random behaviour are not the same. Chaotic systems are deterministic, you can predict the time evolution of the system for a given set of initial conditions. However, even for extremely small changes in the initial conditions for chaotic systems, the time evolution will result in very different trajectories, there being no element of randomness in the trajectory. Chaotic dynamical systems are therefore unpredictable but are not random in nature as any difference in the initial conditions are amplified as the system evolves with time to give enormously different results.

DEFINITIONS

Dynamical system: A system whose behaviour expressed in terms of position, velocity and acceleration may be modelled by a set of differential equations together with a set of initial conditions.

Non-linear dynamical system: A system is non-linear if the set of differential equations contain non-linear terms like, x², xy, e^x, sin(x), etc.

Deterministic system: The solution of the differential equations is unique (one and only one solution) for a given set of initial conditions for deterministic systems.

State: Is a set of specifications of the motion at any time t₀ that is complete enough to determine uniquely the motion of the system at all time t > t₀

Phase space: Phase space is the space of all possible states of a dynamical system.

Fixed point: A dynamical system is in equilibrium at a fixed point where each state variable has a fixed constant value. Fixed points may be classified as stable, unstable or semi-stable.

Orbit (trajectory): The set of points in phase space that satisfy the differential equations that govern the system. For a fixed point, the orbit is a single point.

Periodic orbit: The set of points in the solution of the differential equations repeat themselves and the motion is said to be periodic.

Limit cycle: A closed curve in phase space towards which an orbit may evolve as .

Attractor: The set of points in space towards which the dynamical system evolves as t gets larger. This may be towards a single fixed-point or a limit cycle or an extremely complicated set of points.

Chaotic attractor: If two sets of initial conditions that are almost identical, results in the two orbits diverging exponentially, then the attractor is said to be chaotic. Only non-linear dynamical systems have chaotic attractors. Small differences in the initial conditions lead to widely divergent orbits in phase space making them in practice unpredictable although they are deterministic.

Poincare section: A Poincare map is the intersection of a periodic orbit in the phase (state) space of a continuous dynamical system with a certain lower-dimensional subspace, called the Poincare section. The Poincare section preserves many properties of periodic and quasiperiodic orbits of the original system and has a lower-dimensional state space. So, it is often used for analysing the original system in a simpler way. A Poincare map differs from a recurrence plot in that space, not time, determines when to plot a point.

Bifurcation diagram: A bifurcation is a period-doubling, a change from an N-point attractor to a 2N-point attractor, which occurs when a control parameter r is changed. A Bifurcation diagram is a visual summary of the succession of period-doubling produced as r increases. Bifurcation diagrams are analysed by varying one parameter at a time and keeping others fixed.

SIMPLE PENDULUM (free, damped, driven motion)

We first consider a simple rigid pendulum (zero mass) of length L that is constrained to move along an arc of a circle centred at a pivot point. The angular displacement w.r.t. the vertical is θ, the angular velocity is , and the horizontal displacement is x.

The equation of motion for the pendulum with damping and subjected to a sinusoidal driving force applied at the pivot is

      (1)       

where

           angular displacement [rad]

     t        time [s]

         natural angular frequency [rad.s^-1]

       damping coefficient [s^-1]

      driving strength 

    angular frequency of the driving force [rad.s^-1]

g       acceleration due to gravity  [g = 9.80 m.s^-2]

L      length of pendulum  [m]

To solve equation 1 using the Python function odeint, we need to write this second-order differential equation as the system of two first-order equations

      (2)       

For small amplitude free vibrations of the simple pendulum, its natural period T₀, frequency f₀ of vibration, and angular frequency are

      (3)       

The pivot point is taken as the origin (0,0) and in Cartesian coordinates, to the right is the +X direction and up is the +Y direction. The position of the bob at the end of the pendulum is:

Horizontal displacement at any time t is

      (4)     

Vertical position at the end of the pendulum is

      (5)     

The two fixed-point of the system occur when  and . Hence, the two fix-points of the system are

                                      stable equilibrium point

                    unstable equilibrium point

Any small changes in the damping, driving force strength, or driving force frequency may result in very different trajectories of the pendulum. It is not possible for the motion of the pendulum to be chaotic if the forces acting of the pendulum are only the gravitational force and the damping force.

Chaos can only occur if an external driving force also acts on the system and the driving frequency is less than the natural frequency of oscillation .

If   and as the driving force increases from zero to larger values, you may observe the phenomenon of period doubling. A period doubling bifurcation in a discrete dynamical system is a bifurcation in which a slight change in a parameter value in the system's equations leads to the system switching to a new behaviour with twice the period of the original system. With the doubled period, it takes twice as many iterations as before for the numerical values visited by the system to repeat themselves. A period doubling cascade is a sequence of doublings and further doublings of the repeating period, as the parameter is incremented further and further. Period doubling bifurcations occur in continuous dynamical systems, namely when a new limit cycle emerges from an existing limit cycle, and the period of the new limit cycle is twice that of the old one. However, it is often difficult to find a set of values for the model parameters and initial conditions that demonstrate the phenomenon of period doubling.

SIMULATIONS

Free motion of the pendulum with a small initial displacement   cs_006_01.py

An undamped pendulum with zero driving force acting will vibrate with approximately simple harmonic motion (SHM) when it is given a small disturbance from its equilibrium position. 

Python Console

Initial conditions

  theta(0)/pi = 0.100    omega(0) = 0.0000 rad/s   Damping: b = 0.000

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 0.000   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Results

Peaks t vs x:  Tpeak = 0.671 s   fPeaks = 1.491 Hz

psd:           TPeak = 0.664 s   fPeak = 1.505 Hz

PLOTS

Time evolution plots

Phase space plots

The [2D] plot of  is called the phase space plot or state space plot as the two variables  and  completely defines the state of the pendulum. A phase space orbit is simply the trajectory of the two variables  and  as time evolves.  A closed orbit which is an elliptical attractor always evolves in a clockwise sense. The time to complete one closed orbit is the period of the oscillation. Green dot gives initial state and red dot final state. Lower graph shows the orbit only for the time interval from t­_S to t_F. so that any initial transient vibrations are removed.

Frequency spectrum

The power spectral density psd is plotted against frequency f.

The red vertical line shows the natural frequency f_o and the magenta line for the driving frequency f_D .

The pendulum oscillates at its natural frequency (1.5 Hz) and its motion is described as simple harmonic motion (SHM). The Fourier Transform is calculated by direction integration of the Fourier Transform integral.  For small amplitude oscillation the pendulum corresponds to a linear system.

Free motion of the pendulum with a large initial displacement   cs_006_01.py

The pendulum is released very close to its vertical position (unstable fix-point)  .

Python Console

Initial conditions

  theta(0) = 3.138 rad   omega(0) = 0.0000 rad/s   Damping: b = 0.000

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w = 9.425 rad/s

Driving force

  gamma = 0.000   TD = 1.000 s   fD = 1.000 Hz   wD = 6.283 rad/s

Results

Peaks t vs x:  Tpeak = 1.112 s   fPeaks = 0.899 Hz

psd:           TPeak = 1.087 s   fPeak = 0.920 Hz

Plots

Time evolution plots

Phase space plots

Green dot gives initial state and red dot final state. Lower graph shows the orbit only for the time interval from t­_S to t_F. so that any initial transient vibrations are removed.

Frequency spectrum

The power spectral density psd is plotted against frequency f.

The red vertical line shows the natural frequency f_o and the magenta line for the driving frequency f_D  (not shown as driving strength is zero). Notice the two predominant peaks in the frequency spectrum at 0.90 Hz and 1.50 Hz (natural frequency) and a small peak at 0.31 Hz.

The motion is anything but SHM with a very large initial angular displacement. The period of motion is longer than the natural period or the frequency of vibration is lower than the natural frequency . Notice that the pendulum speeds most of its time near the unstable equilibrium point which occurs at the extreme position of the motion where the angular displacement changes slowly with the angular velocity nearly zero. The time when the angular velocity is zero is the time at which the pendulum changes direction.

It is easy to show if a fixed-point is stable or unstable by considering a small increment away from the equilibrium point, for example,

      à pendulum vibrates around   à unstable

     à pendulum vibrates around   à stable

Free motion of the pendulum with a damping  

Many aspects of damped vibrational motion such overdamped, underdamped, and critical damping can be explored using the code cs_006_01.py.

Python Console

Model Parameters

  theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 0.500

Driving force

  gamma = 0.000   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Time span

  time steps = 5999   tMax = 10.000 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Results

Peaks t vs x:  Tpeak = 0.675 s   fPeaks = 1.482 Hz

psd:           TPeak = 0.632 s   fPeak = 1.581 Hz

Plots

Time evolution

Phase space plot

Frequency spectrum

The motion is attracted to the fixed-point (steady-steady or equilibrium position) of the system .

Forced motions of the pendulum with damping  

The system is excited by some external sinusoidal driving stimulus with amplitude  where   [Python AD] is a strength parameter and f_D is the driving frequency. The response of the system can be studied by changing the input model parameters. Plots can be made for the time evolution of the system, phase space, frequency spectrum, Poincare sections and bifurcation diagrams.

The default parameters are mainly used for each simulation and only the drive strength  is successively increased in the following simulation .

Relatively weak driving strength             The approximation  is valid

Python Console

Model Parameters

  theta(0)/pi = 0.000    omega(0) = 0.0000 rad/s   Damping: b = 3.000

Driving force

  gamma = 0.200   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Time span

  time steps = 5999   tMax = 10.000 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Results

Peaks t vs x:  Tpeak = 0.985 s   fPeaks = 1.015 Hz

psd:           TPeak = 0.945 s   fPeak = 1.059 Hz

Plots

Time evolution

The pendulum motion is to a good approximation SHM.

Phase space

After the initial transient motion, the motion of the pendulum becomes periodic.

Frequency spectrum

The system vibrates at the driving frequency (f_D = 1.00 Hz) and not the natural frequency of vibration (f₀ = 1.50 Hz).

We see that the driver and the response have the same period. Something which intuition from linear problems would say is obvious. The response has two regimes: (1) the decay of an initial transient motion and (2) the steady oscillations at the frequency of the driving signal. The amplitude of the response depends upon the energy balance between the energy supplied by the external driving force and the energy dissipated by the system due to the damping. The phase space plot exhibits a regular orbit which is independent of the initial conditions except for the initial transients which does depend upon the initial conditions. The motion approaches a unique attractor in which the pendulum oscillates sinusoidally with exactly the same frequency as the driving force.

In conclusion for the motion of the linear DDP with a sinusoidal driving force:

(1) There is a unique attractor which the motion approaches, irrespective of the initial conditions applied.

(2) The motion of the attractor is itself sinusoidal with frequency exactly matching the drive frequency.

Weak driving strength             The approximation  is not valid

Python Console

Model Parameters

  theta(0)/pi = 0.000    omega(0) = 0.0000 rad/s   Damping: b = 3.000

Driving force

  gamma = 0.900   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Time span

  time steps = 5999   tMax = 20.000 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Results

Peaks t vs x:  Tpeak = 0.341 s   fPeaks = 2.930 Hz

psd:           TPeak = 0.978 s   fPeak = 1.023 Hz

Plots

Time evolution

After 4 cycles (4 s) the motion settles down to a regular oscillation that looks like sinusoidal with a period equal to the driving frequency. However, the regular oscillations are not sinusoidal since the curve is flatter at the extremes.

Phase space

The motion of the pendulum is periodic with a period equal to the driving frequency (f_D = 1.00 Hz),

Frequency spectrum

Most of the energy supplied by the driving force to the pendulum system excites the fundamental frequency (driving frequency f_D = 1.00 Hz). However, small amount of energy also excites the 3rd harmonic (3.00 Hz) and the 5^th harmonic (5.00 Hz) due to the non-linearity of the DDP system. For the n^th harmonic, the period is given by . Thus, the n^th harmonics will repeat itself n times in one drive cycle. The motion will be periodic at the drive frequency as every harmonic will have cycle back to its original position in one drive period  

So, there is strong evidence (not a proof) that a periodic attractor is approach with a period exactly equal to the driving force.

The boundary between weak and strong driving stimulus is around  for simulations using the default values. For drive strength  the motion of the DDP system becomes very different.

Strong driving strength             The approximation  is not valid

Python Console

Initial conditions

  theta(0)/pi = 0.000    omega(0) = 0.0000 rad/s   Damping: b = 2.356

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.060   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Results

Peaks t vs x:  Tpeak = 0.384 s   fPeaks = 2.607 Hz

psd:           TPeak = 1.009 s   fPeak = 0.991 Hz

Plots

Time evolution

Initially there are dramatic oscillations where the pendulum swings through more than two anticlockwise rotations then swings back to about  radians before it settles down after about 8 drive cycles to sinusoidal like oscillations around 2  radians.  The period of the periodic motion is 1.00 s, which is equal to the period of the driving force. The angular displacement approaches an attractor that oscillates at the frequency of the driving force.

Phase space

The phase space plot shows a single closed orbit after the initial transient period indicating periodic motion of the pendulum centred around . This single closed orbit is called period 1 motion.

Frequency spectrum

The frequency spectrum is characterized by the large peak at the driving frequency (f = 1.00 Hz). The motion of the pendulum becomes complicated because of a number of harmonic oscillations are picked up.

Strong driving strength             The approximation  is not valid

Period 2

Python Console

Initial conditions

  theta(0)/pi = 0.000    omega(0) = 0.0000 rad/s   Damping: b = 2.356

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.073   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Results

Peaks t vs x:  Tpeak = 0.408 s   fPeaks = 2.452 Hz

psd:           TPeak = 1.007 s   fPeak = 0.993 Hz

Plots

Time evolution

There are wild fluctuations in the position of the pendulum during the first 20 s before a regular periodic motion is established. On careful examination of the periodic motion, one will observe that the heights and troughs are not all of the same height, but vary between two distinct heights and this pattern continues indefinitely. So, the strong driving force acting on the system results in a doubling in the period of the pendulum to produce period 2 motion.

Phase space

After the initial wild motion, the phase space plot evolves to two distinct orbits indicating period 2 motion. This behaviour means that the motion no longer repeats itself every drive cycle but every two drive cycles, so the period of the pendulum is twice the drive period. We now have a periodic attractor with period 2.0 s. This phenomenon is known as period doubling

Frequency spectrum (horizontal displacement for periodic motion)

The spectrum is characterized by the two major peaks at frequencies of 1.0 Hz and 2.0 Hz

(f_D = 1.00 Hz). Also, some energy is transferred to the higher harmonics.

The motion of the pendulum is sensitive to small changes in the strength of the driving force. For example, if we increment the drive strength from 1.073 to 1.075 the attractor shifts from an orbit around -2  to around +4  with period 2 motion.

Although the attractor has period 2, the dominant behaviour is still clearly period 1.

Strong driving strength             The approximation  is not valid

Period 3

There is now an attractor for which a subharmonic term is dominant and the motion settles down to an attractor that repeats itself every 3 drive cycles and hence period 3 (T = 3.0 s).

Initial conditions

  theta(0)/pi = 0.000    omega(0) = 0.0000 rad/s   Damping: b = 2.356

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.077   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

Results

Peaks t vs x:  Tpeak = 0.590 s   fPeaks = 1.695 Hz

psd:           TPeak = 0.984 s   fPeak = 1.017 Hz

Time evolution

The motion quickly settles down to period 3 motion (T = 3.00 s)

Phase space

A single orbit is established after the initial transient motion.

Frequency spectrum (angular velocity )

The subharmonic term becomes the dominant term in the frequency spectrum .

Initial conditions

For a linear oscillator, with a given set of parameters there is a unique attractor that is independent of the initial conditions, the eventual motion will be the same once the transients have decayed away.  This is not the case for a non-linear system as shown below. The initial condition for the blue curve is , whereas the initial condition for the red curve is  .

For the motion with the initial condition  the period of oscillation is 3 s. When the initial condition was   the period is actually 2 s with the peaks and troughs having slightly different heights

Period doubling cascade

Simulation parameters:

As the driving stength  is incremented, the motion changes to different attractors as you get a sequence of bifurcations leading to a period-doubling cascade:

         = 1.0600, T = 1.00 s, pattern repeats itself every one drive cycles

         = 1.0780, T = 2.00 s, pattern repeats itself every one drive cycles

         = 1.0810, T = 4.00 s, pattern repeats itself every one drive cycles

         = 1.0826, T = 8.00 s, pattern repeats itself every one drive cycles

         = 1.0828, T = 16.00 s, pattern repeats itself every one drive cycles

   Motion settles down to steady oscillations at the same frequency of the drive excitation. The attractor has a period of 1.00 s. The expanded view shows that all the peaks have the same height.

  Motion settles down to steady oscillations at the same frequency of the drive excitation. But now, the peaks are not all the same height. The maxima alternate between two fixed heights, so the attractor has a period of 2.0 s.

   Motion settles down to steady oscillations at the same frequency of the drive excitation. But now, the peaks are not all the same height. The maxima alternate between four fixed heights, so the attractor has a period of 4.0 s.

   Motion settles down to steady oscillations at the same frequency of the drive excitation. But now, the peaks are not all the same height. The maxima alternate between eight fixed heights, so the attractor has a period of 8.0 s. the pattern repeats itself every 8 drive cycles.

   Motion settles down to steady oscillations at the same frequency of the drive excitation. But now, the peaks are not all the same height. The maxima alternate between sixteen fixed heights, so the attractor has a period of 16.0 s. the pattern repeats itself every 16 drive cycles.

It is remarkble that the phenomenum of period-doubling is observed in many non-linear physical systems as a control parameters is incremented. In all such systems, the period-doubling cascade occurs in the same way, a circumstance known as universality.

The Feigenbaum Number and Universality

The period doubling occurs more and more frequently as  is increased. To examine this in more detail, we need to look at threshold values (bifurcation values) of   for period doubling to occur. Finding the bifurcation values is diffiecult as you may need at least 5 significant figures. For our intial values, the first four bifurcation values are:

                     1 à 2    

                     2 à 4    

                    4 à 8     

                    8 à 16   

We can define the Feigenbaum number as

For many non-linear systems, the Feigenbaum number is universal and, in the limit, as

In our example of the DDP system

If we continue this sequence of period doubling, we will reach a limit, to give a critical value for

For  the period-doubling terminates and chaos sets in. So, the period-doubling cascade is often called the route to chaos (note: there are other routes to chaos without the period-doubling occurring).

CHAOTIC MOTION

For drive strengths greater than about a critical value of , the solution is not even periodic at all! In linear theory, a damped system driven periodically must eventually respond periodically at the driving frequency and the oscillations independent upon the initial conditions.  In this case for drive strengths greater than the critical value, the DDP nonlinear system will oscillate forever without ever repeating – it is chaotic. The motion of the pendulum becomes erratic and unpredictable since very small changes in any parameters such as time step, initial values may result in very different trajectories.

The case for   with initial conditions , the DDP is obviously trying to oscillate at the drive frequency f_D = 1.00 Hz. Nevertheless, the actual oscillations wander around erratically without ever repeating themselves. This erratic and non-periodic motion is one of the chief features of chaotic motion.

CHAOS and sensitivity to initial conditions

The other defining feature of chaos is that a trajectory is extremely sensitive to the initial conditions. For our DDP system with  two orbits with slightly different initial conditions are almost identical, then they follow different trajectories. This makes it an impossibility to predict the trajectories for chaotic motion.

Lyapunov Exponents    cs_006_02.py

Consider two identical pendulum which are released with different initial conditions for the angular displacement . To a reasonable approximation, the difference in the angular displacements  will vary as

where the parameter  is known as the Lyapunov exponent. and k and k’ are constants.  when plotted on a log scale against time will be a straight line with the slope of the line equal to the value of the Lyapunov exponent.  If , then  and the motions will converge exponentially.  If , then  and the motions will diverge exponentially. 

Consider two motions of the pendulum with a small driving force   with different initial conditions  . The upper graph shows the angular displacements as functions of time and the lower graphs, the difference in the angular displacements   where  is plotted on a linear and a logarithmic scale.

It is clear from log scale plot, the maxima in the log scale for  decrease linearly, hence  decays exponentially, dropping 6 orders of magnitude in the first six drive cycles . In the linear regime, the separation  of two identical DDPs with different initial conditions, decreases exponentially with time . Linear oscillators are insensitive to the initial conditions. So, to make accurate predictions, you only need to know the initial conditions to the same accuracy.

With a moderate drive strength   the motions are periodic and their trajectories converge with  decaying exponentially. For the motion governed by a moderate driving signal, the sharp dips in the   plot occur when one of the pendulums reaches a turning point, following which  will vanish since the two pendulums will cross each other. You will notice that  decreases steady with time such that  (ignoring the dips). This means that the motion of the pendulum is predictable. Knowledge of the motion of the first pendulum enables you to predict the motion of the second pendulum even though you don’t know the second pendulum’s initial conditions when the forcing is small to moderate.

Initial condtions:  

Two identical DDPs with a strong driving force   that were released with an extremely small difference in initial angular displacements. The motion for the large excitation initially shows that the two motions are almost identical then suddenly they diverge as  increases rapidly and irregularly with time . Hence, the motions of the two pendulums diverge exponetially from one another with time. So, for chaotic damped driven motion of the pendulum, we simple do not know where the pendulum will be in the future, since even for minute differences in the initial conditions it will lead to trajectories that are very different.

For about 12 s, the motion of both pendulums is erratic but follow similar orbits and then the orbits diverge exponentially then levels out.

The above results for pendulums which start with nearly identical initial conditions indicate that for small or medium forcing, the motions will converge exponentially whereas for high forcing, the trajectories diverge exponentially. The pendulum system while obeying deterministic laws may still exhibit irregular and unpredictable behaviour due to an extreme sensitivity to initial conditions.      

BIRFURCATION DIAGRAM     cs_006_03.py

The purpose of a bifurcation diagram is to show in a single plot, the changing periods, alternating periodicity, and chaos as the drive strength   varies. It is a plot of  or a plot of .

The required steps in plotting a bifurcation diagram are:

1.    Choose a large number of evenly spaced values for    from   to  .

2.    Choose the initial conditions   and  .

3.    Solve the equation of motion for DDP for each value of  for t = 0 to t = t_max.

4.    Check for periodicity or non-periodicity when all the transience has disappeared by examining  or  from t_min to t_max in one-cycle intervals with period T.

5.    Plot these values of   or  from t_min to t_max against each value of .  For larger values of  the pendulum can start to make many revolutions, so it is necessary to restrict   to the range  .  It is often better to plot rather than  so that you don’t have to worry about the multiple rotations of the pendulum.

The bifurcation diagram is very dependent upon choice of parameters and initial conditions. The following diagrams show the period doubling cascade effect for initial conditions . It took about 5 to 10 minutes to produce each bifurcation diagram using Python (Sypder).

Bifurcation:

PHASE SPACE ORBITS: POINCARE SECTION    cs_006_01.py   cs_006_04py

A Poincare section is another way to visualise the motion of chaotic systems. The Poincare section is a simplified view of a phase space orbit. The phase space or state-space is the plot of . As time passes, the point  traces out an orbit in the phase-space plot. Such plots give a very clear picture of the motion of the system.

cs_006_01.py

Initial conditions

  theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 2.356      t_Max = 10 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 0.600   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

The orbit spirals inwards in a clockwise direction and rapidly approaches the period 1 attractor and then continually cycles around the attractor. 

If the initial angular displacement is  then the orbit spirals outward to the elliptical attractor. After the transient motion has decay away, the orbit follows the elliptical attractor in a clockwise direction.

So, for large t values, the two orbits become identical and one can not distinguish between the two different initial conditions. The phase-space plot gives a clearer picture of the approach to the attractor than the time evolution plots.

cs_006_01.py

Initial conditions

  theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 2.356      t_Max = 40 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.078   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

After the transient motion has decayed way, the orbit is composed of two distinct loops, and the attractor is period 2.

cs_006_01.py

Initial conditions

  theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 2.356      t_Max = 40 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.081   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

After the transient motion has decayed way, the orbit is composed of four distinct loops, and the attractor is period 4.

cs_006_01.py    CHAOS

Initial conditions

  theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 2.356      t_Max = 40 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.105   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

There is no closed orbit, the motion never repeats itself, so the motion is chaotic.

cs_006_01.py    ROLLING MOTION

Initial conditions

  theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 2.356      t_Max = 6 s

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

  gamma = 1.400   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

For large values of the drive strength, we can a rolling motion where the pendulum swings through a complete revolution each drive cycle. After the initial transient motion, the motion of the pendulum is perfectly periodic as the motion in each loop is identical.

POINCARE SECTION             cs_006_04.py

There is a way to study chaotic motion that is better than simply plotting the trajectory in phase space because after many cycles it contains too much information to be useful. Consider the phase space plot where points are not plotted at every time step but only at times given by

Such a phase space plot is called a Poincare section. If the pendulum oscillates at the driving frequency, then only one point will appear in the Poincare section. If the oscillation has twice the frequency of the driving force, then the Poincare section will have two points. If the motion is not periodic and chaotic the Poincare section will consist of a pattern of points called the attractor. The attractor has a structure that is frequently beautiful even though the motion is unpredictable and chaotic, yet at the same time preserve a coherent global structure. So, in the Poincare section plot, the orbit is not drawn but only points at one drive cycle interval. For periodic motion the Poincare section is not useful, but it is very useful for visualisation of chaotic motion.

Period 1 oscillations                 cs_006_04.py

For   the motion is periodic with the period equal to the drive period (1.00 s).    

The phase space plot (red: orbit after transient motion decayed away) and Poincare section for  . Because the oscillation is described as period 1 motion, there is only one dot in the Poincare section.

Period 2 oscillations     

For   the motion is periodic with the period equal to twice drive period (1.00 s).

The red trajectory shows the orbit after the transient motion has decayed away. The motion is period 2 since there are two loops that make up one period of the oscillation and two dots are show on the Poincare section.

Poincare section for chaotic motion

Initial conditions     

      theta(0)/pi = -0.500    omega(0) = 0.0000 rad/s   Damping: b = 1.178   

Free vibration

               T0 = 0.667 s   f0 = 1.500 Hz    w0 = 9.425 rad/s

Driving force

gamma  = 1.500   TD = 1.000 s   fD = 1.000 Hz    wD = 6.283 rad/s

The pendulum undergoes an erratic rolling motion making many complete revolutions in one direction and then in the other direction, but never repeating itself.  The phase space plot is not useful because of the entanglement of the orbit as the pendulum swings through complete revolutions in one direction, then the other.

The Poincare section and enlargements

The Poincare section contains a subset of all the points of the phase space orbit, and it is impossible to know what this subset of points will look like, but we are able to compute and display it. The Poincare section often gives a very elegant picture by plotting this subset of points which are one drive cycle apart. The Poincare section is simply not a figure but a fractal where one finds further structure by enlarging sections of the Poincare section.  For example, zooming in a “tongue” is actually made up of many tongues. If you had plotted enough points, you could keep zooming in a finding a repetition of the tongue structure. The Poincare section took about 2 hours to compute and plot. This self-similarity is one of the features of a fractal. for a chaotic system which can be visualised as a fractal, the long-term motion is said to be a strange attractor.