Chaos

Chaos#

Author:
Date:
Time spent on this assignment:

Note: You must answer things inside the answer tags as well as questions which have an A:.

In this project we are going to learn about pendulum and chaos. We will start with a single pendulum, examining phase space diagrams and animating trajectories, and moving toward a double pendulum.

Exercise 1. Single Pendulum#

List of collaborators:
References you used in developing your code:

You will use the following animation code for exercise 1

Help with animateMe: There are two animate functions that accept positions (the angles, technically) and l1 which is the length of the pendulum. The positions argument needs to be passed as a nested list [positions] and the positions array needs to be 2-dimensional. You can check the shape of your array using positions.shape to ensure it is in the right format of (N,1), which is what the animation functions expect. If your positions array is 1D, you can reshape it using positions[:,np.newaxis].

a. A Single pendulum#

When working with pendula, instead of keeping track of the position \((x,y)\), we instead are going to keep track of the angle \(\theta\).

In this exercise, we are going to simulate a pendulum with a rigid rod \(L\) and a fixed mass \(m\) at the end of it. We can define the position of the pendulum by the angle it makes with respect to hanging down (at \(\theta=0\)).

The “force” (angular acceleration) of a single pendulum is

\[ \alpha(\theta) \equiv \frac{d^2 \theta}{dt^2}= -B\sin(\theta), \]

where \(B\) (on Earth) is the gravitational constant \(g/L\).

🦉Write a function that uses the midpoint method to simulate a pendulum:

def simulate_single_pendulum(p, s):
    """
    Simulate a single pendulum 
    p: SinglePendulumParams that defines the system
    s: SinglePendulumSimulation that defines the particular simulation parameters

    Returns: 
       times: np.ndarray (nstep,)
       thetas: np.ndaray (nstep,)
       omegas: np.ndarray (nstep,)
    """

NOTE: If you have had trouble getting the midpoint method to work, please ask for help!

We first want to run a simulation with the following parameters:

B = 9.8 (this is the force of gravity)
L = 1.0
m = 1.0 (actually m doesn’t matter until we add driving forces later)
A = C = omega_drive = 0 (these will be driving/damping which we will treat later) You can make a parameter object as follows:

p1 = SinglePendulumParams(L = 1, B = 9.8)

Run a simulation with no initial velocity and an initial angle of \(\theta=0.6\) and a \(dt=0.01\) for a time \(T=10\). You can set this up by

s1 = SinglePendulumSimulation(dt=0.01, T=10, theta_0=0.6, omega_0=0.0)

You then should be able to generate the trajectory using your function by calling

times, thetas, omegas = simulate_single_pendulum(p1,s1)

Plot

\(\theta\) vs. \(t\)
\(y\) vs. \(x\) (include the origin on this plot, and equalize the axes scale)
a phase space diagram (\(\omega\) vs \(\theta\))
and animate the pendulum.

A phase space diagram is a graph displaying the position and velocity of an object on the abscissa (x axis) and ordinate (y axis). Since an undamped, undriven harmonic oscillator moves with

\[\theta(t) = A\cos(\Omega t + \phi)\hspace{2cm} \omega(t)=-A\Omega\sin(\Omega t+\phi)\]

its phase trajectory will be an ellipse with axes of lengths \(2A\) and \(2A\Omega\) (unrelated to the \(A\) in Params). Confirm that you get this.

from dataclasses import dataclass
@dataclass
class SinglePendulumParams:
    """
    Define the parameters of the single pendulum system
    """
    L: float   # length of the pendulum
    A: float = 0.0 # damping
    B: float = 9.8 # restoring force
    C: float = 0.0 # driving force
    omega_drive: float = 0.0 # driving frequency 

@dataclass
class SinglePendulumSimulation:
    """
    Define the simulation parameters for the single pendulum system
    """
    dt: float = 0.01 # time step
    T: float = 20.0 # total time of the simulation
    omega_0: float = 0.0 # initial angular velocity
    theta_0: float = 0.0 # initial angle

Q: Does the phase space diagram match the analytic expectations? Explain how you know.

A:

b. Pendula and different starting positions#

Recall in your intro mechanics class (Physics 211/etc) that you know the analytic solution for a pendulum is

\[ \theta(t) = \theta_0 \cos(\Omega t) \]

However, that solution only works for small angles.

🦉Let’s quickly check this, by looking at the between the small-angle analytic formula and the numerical solution for \(\theta_0 \in \{0.1,0.3,1\}\).

For each value of \(\theta_0\), plot the phase space diagram using the intro mechanics solution and the numerical solution on top of one another.

It may help to use subplots; for example, you can do:

fig, axes = plt.subplots(1,3, figsize=(9,3)
axes[0].plot(thetas1, omegas1, label='numeric')
axes[0].plot(thetas1_analytic, omegas1_analytic, label='analytic')
axes[0].legend()
axes[0].set_title("$theta_0$=0.1$)
axes[1].plot(thetas2,omegas2)

etc. At the end of it if your plots are too closely spaced, you can use plt.tight_layout()

###ANSWER HERE

Q: Describe the difference between the small-angle formula and numerical result as a function of the amplitude.

A:

Q: Which solution do you expect to be more accurate? Can you explain the trend from above?

A:

c. Damped and Driven and phase plots#

The behavior of a “simple” pendulum becomes less simple when you add damping and drive the system.

Suppose we make a pendulum from an object of mass \(m\) at the end of a rigid, massless rod of length \(L\). The rod is suspended from a bearing which exerts an angular velocity-dependent drag force on the pendulum. Because we are using a rod, not a string, the pendulum can swing in a full circle if it is moving sufficiently quickly. A motor attached to the pendulum provides a sinusoidally varying torque. Gravity acts on the pendulum in the usual way.

The equation of motion for the pendulum is then

\[\begin{split} \begin{align} \tau = I\alpha &= mL^2 \frac{d^2\theta}{dt^2}\\ &=-mgL\sin\theta - \beta \frac{d\theta}{dt} + \gamma \sin \Omega_\mathrm{ext} t \end{align} \end{split}\]

The first term to the right of the equal sign is the torque due to gravity. The second term is the velocity-dependent damping, while the third term is the external driving torque. The symbols \(\beta\) and \(\gamma\) represent constants.

Do some algebra:

\[\frac{d^2\theta}{dt^2} = -\frac{g}{L} \sin\theta -\frac{\beta}{mL^2}\frac{d\theta}{dt} + \frac{\gamma}{mL^2}\sin\Omega_\mathrm{ext} t\]

or

\[\frac{d^2\theta}{dt^2}+A\frac{d\theta}{dt} + B\sin\theta = C\sin\Omega_\mathrm{ext} t\]

where \(A\) and \(B\) are positive constants. We can rewrite this as a pair of first order equations:

\[\frac{d\theta}{dt}=\omega \hspace{3cm} \frac{d\omega}{dt} = -A\omega -B\sin\theta+C\sin\Omega_\mathrm{ext} t\]

🦉Update simulate_single_pendulum to calculate (and graph) \(\theta\) vs. \(t\) and \(\omega\) vs. \(t\) for the pendulum. You may find it convenient to write your expression for the force in this way:

def force_single_pendulum(t,theta,omega,p):
    """
    t : current time (s)
    theta: current angle
    omega: current angular velocity
    p: SinglePendulumParams 
    Returns:
    d omega / dt 
    """
    F = -p.A*omega - p.B*np.sin(theta) + p.C*np.sin(p.omega_drive * t) 
    return F

Make sure when you make your plots that all values of \(\theta\) are \(-\pi \leq \theta \leq \pi\). When you produce your graphs make sure that you wrap things so that those are the only values that you see on your graph. You may want to use plt.scatter so that there are not annoying lines.

Run your program with the following parameters and be sure to include all the plots in your Jupyter notebook.

Parameter	Version 1	Version 2	Version 3
\(A\) (damping)	0	0.1	0.1
\(B\) (restoring force)	0.1	1	1
\(C\) (external driving force)	0	0.1	2
\(m\)	1	1	1
\(\Omega_\mathrm{ext}\) (driving frequency)	N/A	2	1.2
\(\theta_0\)	0	0	0
\(\omega_0\)	0.1	0.1	0.1
\(t_{max}\)	120	120	500
\(dt\)	.01	.01	.01

Version 1 is an undamped, undriven pendulum. The amplitude of motion is small enough so that it behaves very much like a harmonic oscillator with oscillation frequency \(\Omega=\sqrt{0.1}\). Version 2 is lightly damped, and driven at twice its natural frequency. As the initial motion (which has \(\Omega=1\)) dies away, the driving force will cause it to oscillate at the same frequency as the driving frequency.

It takes longer for the pendulum to settle down than it did in case 2.

If we had tuned the parameters just right, we might have found a set for which the motion was chaotic, never settling into a periodic motion.

### ANSWER HERE

d. Simple pendulum phase space trajectory#

The phase trajectories of driven nonlinear oscillators can be much more complicated than this. 🦉Please generate phase space diagrams for the simple pendulum using the parameters for version 2 and version 4 (shown below again). Note that \(t_\textrm{max}=500\)

Parameter	Version 2	Version 4
\(A\) (damping)	0.1	0.1
\(B\) (restoring force)	1	2
\(C\) (external driving force)	0.1	2
\(m\)	1	1
\(\Omega_\mathrm{ext}\) (driving frequency)	2	1.2
\(\theta_0\)	0	0
\(\omega_0\)	0.1	0.1
\(t_{max}\)	500	500
\(dt\)	.01	.01

###ANSWER HERE

If you want, go ahead and animate these pendula

Exercise 2. Double Pendulum#

List of collaborators:
References you used in developing your code:

Single (rigid) pendula are somewhat boring. We can compute most things about them analytically. A double pendulum on the other hand has a much richer behavior. In this exercise you are going to simulate a double pendulum.

Show code cell content

Hide code cell content

def animate_double_pendulum(positions, l1, l2): 
    """
    Animate a trajectory
    positions: list of np.array, each of shape (T,2) where T is the number of time steps and 2 corresponds to the two angles of the pendulum
    l1: length of the first pendulum
    l2: length of the second pendulum
    returns: a matplotlib animation object that can be displayed in a Jupyter notebook by running HTML(ani.to_jshtml())
    """
    positionArray=[]
    for position in positions:
        theta1=position[:,0]
        theta2=position[:,1]
        x1=l1*np.sin(theta1)
        y1=-l1*np.cos(theta1)
        x2=l1*np.sin(theta1)+l2*np.sin(theta2)
        y2=-l1*np.cos(theta1)-l2*np.cos(theta2)
        l=len(x1)
        positions=np.zeros((l,4))
        positions[:,0]=x1
        positions[:,1]=y1
        positions[:,2]=x2
        positions[:,3]=y2
        positionArray.append(positions)
    fig, ax = plt.subplots()
    x_min=np.min([np.min(list(positions[:,0])+list(positions[:,2])) for positions in positionArray])*1.1
    x_max=np.max([np.max(list(positions[:,0])+list(positions[:,2])) for positions in positionArray])*1.1    
    y_min=np.min([np.min(list(positions[:,1])+list(positions[:,3])) for positions in positionArray])*1.1
    y_max=np.max([np.max(list(positions[:,1])+list(positions[:,3])) for positions in positionArray])*1.1
    y_max=np.max([y_max,0.1])
    y_max=np.max([y_max,x_max])
    x_max=y_max    
    y_min=np.min([y_min,x_min])
    x_min=y_min
    ax = plt.axes(xlim=(x_min, x_max), ylim=(y_min, y_max));
    ax.plot(['0'],['0'],'o')
    lines=[]
    colorWheel=['g','b','r']
    for idx,positions in enumerate(positionArray):
        lines.append((ax.plot([], [],'o', color = colorWheel[idx%len(colorWheel)]))[0])
        lines.append((ax.plot([], [],'-', color = colorWheel[idx%len(colorWheel)]))[0])
        lines.append((ax.plot([], [],'o', color = colorWheel[idx%len(colorWheel)]))[0])
        lines.append((ax.plot([], [],'-', color = colorWheel[idx%len(colorWheel)]))[0])
        
        
    def update(i, positionArray,lines):
        for idx,positions in enumerate(positionArray):
            lines[4*idx+0].set_data([positions[i,0]],[positions[i,1]])
            lines[4*idx+1].set_data([0,positions[i,0]],[0,positions[i,1]])
            lines[4*idx+2].set_data([positions[i,2]],[positions[i,3]])
            lines[4*idx+3].set_data([positions[i,0],positions[i,2]],[positions[i,1],positions[i,3]])
        return lines
    ll=1
    ani = animation.FuncAnimation(fig, update, len(positionArray[0]), fargs=[positionArray, lines],
                      interval=20, blit=True,repeat=False)
    plt.close()
    return ani

@dataclass
class DoublePendulumParams:
    """
    Define the parameters of the double pendulum system
    """
    L1: float   # length of the first pendulum
    L2: float   # length of the second pendulum

    m1: float = 1.0 
    m2: float = 1.0 
    g: float = 9.8 # restoring force

@dataclass
class DoublePendulumSimulation:
    """
    Define the simulation parameters for the double pendulum system
    """
    omega_0: np.ndarray # initial angular velocity (2,)
    theta_0: np.ndarray # initial angle (2,)

    dt: float = 0.01 # time step
    T: float = 20.0 # total time of the simulation

a. Simulating double pendula#

Modify your code to simulate a double pendulum. Now instead of representing the location of your pendulum with one \(\theta\) you need to use two \(\theta\). (Be careful you are dealing with the mass correctly). Again, this can be mainly accomplished by playing with your initial conditions.

Now, you need to get the force to work. The force for a double pendulum is somewhat complicated. To compute the force you either need to:

derive the force carefully
- This is a bit annoying because part of the force involves the tension in the rod (i.e. the first mass has a force of \((-T_1 \sin(\theta_1) +T_2 \sin(\theta_2), T_1 \cos \theta_1 - T_2 \cos \theta_2 - m_1 g) ) \) Because we don’t know what \(T\) is, we need to do some algebraic manipulation to be able to write \(\frac{d^2\theta}{dt^2}\) as a function of \(\omega=\frac{d\theta}{dt}\) and \(\theta\) (which is what we need)
use the Euler-Lagrange equations. This ia a better approach for systems with constraints and is something you will learn in classical mechanics.

At the moment we will just give you the effective force¹:

def force_two_pendulum(t1, t2, omega1, omega2, p):
    """"
    t1: angle of the upper pendulum
    t2: angle of the lower pendulum
    omega1: angular velocity of upper pendulum
    omega2: angular velocity of the upper pendulum
    p: a DoublePendulumParams object
    Returns:
    omega1_dot, omega2_dot : angular acceleration
    """
    f1=-p.l2/p.l1 *(p.m2/(p.m1+p.m2))*omega2**2*np.sin(t1-t2)-p.g/p.l1*np.sin(t1)
    f2=p.l1/p.l2 * omega1**2 * np.sin(t1-t2)-p.g/p.l2*np.sin(t2)
    alpha1=p.l2/p.l1*(p.m2/(p.m1+p.m2))*np.cos(t1-t2)
    alpha2=p.l1/p.l2*np.cos(t1-t2)
    den=(1-alpha1*alpha2)
    omega1_dot = (f1-alpha1*f2)/den
    omega2_dot = (-alpha2*f1+f2)/den
    return omega1_dot, omega2_dot

Note that omega1/omega2 here mean the velocity or \({d\theta_i}/{dt}\), \(i=1,2\) respectively. Also, in the notation above \(t1=\theta_1\) and \(t2=\theta_2\), not to be confused with the tensions!

Put this force in your code.

🦉Run with the following parameters:

    m1:1.0,
    m2:1.0,
    L1:1.0,
    L2:1.0,
    g:9.8,

    T :15.0,
    dt:0.01,
    theta_0 : np.array([1,1+0.11]),
    omega_0 : np.array([0.0,0.0])
     

and make a phase space plot of the system for both the first mass and the second mass (on the same plot). Also go ahead and animate it.

Then switch to use the initial position of theta_0 = np.array([3.14,3.14+0.11]),. How do these compare?

In the animation, you may not want to use all the points.

### ANSWER HERE

Q: Does the inner bob (\(m_1\)) act like a simple pendulum? Use the phase space diagram in your answer

A:

Q: How does the behavior of the pendulum change compared to a simple pendulum as the initial angle is increase?

A:

b. Chaos#

Our goal now is to see that this double pendulum is chaotic. Chaos means that if the system starts with very similar initial conditions, that the result will look very different in the future (long times). Go ahead and run a series of 10 initial conditions except that each one differs from the previous one by an initial angle on the second pendulum of \(10^{-4}\), Animate them all simultaneously to see what happens.

Hint: Your code should have something like this

p = DoublePendulumParams(m1=1.0, m2=1.0, L1= 1.0, L2=1.0, g=9.8)
for i in range(10):
    theta2 = 3.14+0.11+i*1e-4
    s = DoublePendulumSimulation(T=15.0, 
                                 dt=0.01, 
                                 theta_0 = np.array([3.14, theta2]),
                                 omega_0 = np.array([0.0,0.0]) 
                                 )

    ts, thetas, omegas = simulate_double_pendulum(p,s)
    # collect and plot data

###ANSWER HERE

Acknowledgements:

Ex 1: George Gollin (original); Bryan Clark, Ryan Levy (modifications)
Ex 2: Bryan Clark (original)