Chaos in the
Pendulum

We have already been introduced to some of the features of the simple pendulum. In the laboratory we will investigate the behavior of the simple pendulum as we make our model of the pendulum more realistic by including friction and a driving force.

1. Making a Realistic Model.

The equation describing the  physical pendulum can be written in terms of torques as

                                                  τ = Iα =  (Subscript[F, p]) r
                                                
where (Subscript[F, p]) is the perpendicular force acting on the pendulum, r is the moment arm, I is the moment of inertia, and α is the angular acceleration (d^2)θ/(dt^2). For a uniform rod of length L pivoted at one end and swinging under the influence of gravity only one obtains the following expression

                                                I (((d^2) θ)/(dt^2)) = (-mg sinθ ) (Subscript[L, cm]) 

where (Subscript[L, cm]) is the length from the pivot on the end to the center of mass (CM) of the pendulum, m is the mass, and  θ is the angle with respect to the vertical. To make our model more realistic we will include torques caused by friction and the existence of an outside 'driving' torque. This last torque could be caused by a variety of sources (e.g., gusts of wind, a varying electric field acting on a charged pendulum) and will open up a rich array of interesting behavior in the system.

The torque created by friction typically depends on the speed of the object through a viscous medium. In the appropriate velocity range one finds this torque can be treated as a force acting on the center of mass and can be written as

                                         (Subscript[τ, f]) = -(Dd(θ/dt)) (Subscript[L, cm])
                                         
where D is the drag coefficient detemined by the nature of the physical pendulum and (Subscript[L, cm]) is the pivot-CM distance of the pendulum. The minus sign is critical since this torque must oppose the motion.

We now add a time-varying 'driving' term which corresponds to some outside influence acting on the system. This torque has the form

                                       (Subscript[τ, D]) = (Subscript[L, cm])(Subscript[F, D])sin(Ωt)
                                       
 where  (Subscript[F, D]) is the amplitude of the force and Ω is the angular frequency of the driving force.
 
 a. Starting with the pieces above find the full differential equation that will describe the physical pendulum moving uncer the influence of gravity, friction, and a time-varying driving force.
 
 b. Using centered formulas for any derivatives in your expression find a recursion relationship that will enable you to numerically solve your differential equation. What information do you need to get the algorithm started? Clearly show how you would generate the necessary information to get the algorithm started.

2. Solving the Equation

Write a code that will solve the model of the physical pendulum you developed in Section 1. The pendulum is a uniform cylinder that is pivoted at one end. Be sure to account for the fact that -π<θ<π. The following inital conditions and parameters apply.

\(TraditionalForm\`t\_0\  = \ \(0.0\ s\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ θ\_0\  = \ \(25\^O\ \ \ \ \ \ \ \ \ \ \ \ \ \ ω\_0\  = \ 0.0\ rad/s\)\)\)

Length of pendulum:  3.0 m (it's big)
Drag coefficient: 0.5 kg-m/s
Mass of pendulum: 30 kg
Driving force: 250 N
Angular frequency of driving force: 0.8 rad/s

Integration step size: 0.002 s

3. Testing the Code

Test your code by generating a plot of angular position θ with respect to time. If you set the driving force and the drag coefficient to zero you should recover the behavior of the simple physical pendulum. Check that the amplitude you observe is consistent with the initial conditions for a simple physical pendulum. Next, turn on the drag coefficient. What do expect the pendulum to do? Does your code reproduce this behavior?

4. Investigating the Physics

Once you have a working code consider the following questions.

a. How does the time dependence of θ in your new model differ from what you expect for the simple pendulum? 

b. The most salient feature of dynamical chaos is sensitivity to initial conditions. This means that the same system with just slightly different initial conditions will follow a dramatically different trajectory. In fact, 'adjacent' sets of initial conditions give rise to exponentially diverging trajectories. Try generating a plot of angular position versus time for (Subscript[θ, 0])  =  (26^o). Does your model exhibit this feature? Explain.

c. Modify your code to produce a phase-space portrait of the trajectory. What do expect this figure to look like for the simple physical pendulum? How does your plot differ? What information could you extract from the phase-space portrait of the pendulum?

d. A feature of the plots above is that the pendulum can only swing through a range of angles of -π to π. Consider the example below that demonstrates how one could change this feature. Employ this method for your plot of angular position versus time. Do your results make sense?

5. Some Necessary Commands.

If the pendulum swings outside the range of angles of -π to π one can use the 'If' function in Mathematica to correct for this feature. The code fragment below demonstrates the use of this method. Compare the plots with and without the restriction on the angular range.

\(\(\( (*\ 
      constants\ *) \)\(
\)\(
\)\(\(tmin\  = \ 0.0;\)

    \(tmax\  = \ 5.0;\)

    \(tstep\  = \ 0.1;\)

 (*\ 
      a\ function\ for\ \(\(theta\)\(.\)\)\ *) 


    \(angle[t_]\  := \ t^2\  - \ 4;\)

 (*\ 
      Make\ a\ table\ in\ which\ θ\ only\ covers\  - π\ to\ \(\(π\)\(.\)\)\ *) 


    \(t1\  = \ 
        Table[{t, 
	
            theta\  = \ angle[t]; 
	
            If\ [angle[t]\ \  > \ π, \ 
              theta = angle[t]\  - \ 2*π]; 
	
            If\ [angle[t]\ \  < \ π, \ 
              theta = angle[t]\  + \ 2*π]; 
	
            theta
}, 
{t, tmin, tmax, 
            tstep}
];\)

 (*\ 
      plot\ \(\(it\)\(.\)\)\ *) 


    ListPlot[t1, 
	\ PlotJoined -> True, Frame -> True, 
	
      FrameLabel -> {"t", "Angle", "", ""}]\)\)\)

To make a 3-dimensional plot of (t,θ,ω) you will some new commands. These are not in the default version of Mathematica so they have to loaded in. This is done below.

<< 
      Graphics`Graphics3D`

\(\(\( (*\ Make\ a\ 3  D\ table\ of\ time, \ angular\ position, \ 
      and\ angular\ velocity\ and\ plot\ \(\(it\)\(.\)\)\ *) \)\(
\)\(
\)\(\(t2\  = \ 
        Table[{t, Cos[2*t], \(-2.0\)*Sin[2.0*t]}, {t, 0, 10, 
            0.1}];\)

    ScatterPlot3D[t2, 
	
      PlotJoined -> True, 
	
      PlotStyle -> RGBColor[1, 0, 0], 
	
      PlotRange -> {{0, 10}, {\(-1\), 1}, {\(-2\), 
            2}}, 
	\ \ BoxRatios -> {1.5, 1, 1}, 
	Axes -> True, 
	
      AxesLabel -> {"t(s)", "θ(rad)", "ω(rad/s)"}, 
	PlotLabel -> "Scatter Plot", 
	
      ViewPoint -> {2.002, \ \(-2.494\), \ 1.104}
]\)\)\)