Numerical project 4

This 25 point project is due next Wednesday, April 20. Your mission is to generate an animation that looks something like so: I'll get you started. The individual frames can be generated with the following code.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Polygon

def pic(theta):
    x_end,y_end = (np.cos(-np.pi/2+theta), np.sin(-np.pi/2+theta))
    plt.plot([0,x_end],[0,y_end], 'k', linewidth=5)
    plt.plot(x_end,y_end,'ro', markersize=10)
    plt.plot([0,0],[0,-1], 'k', linewidth = 5, alpha=0.1)
    plt.plot(0,-1,'go', markersize=10, alpha=0.6)
    ax = plt.gca()
    ax.add_patch(Polygon([[-0.3,0],[0.3,0],[0.3,0.1],[-0.3,0.1]], facecolor='k'))

    ax.set_xlim(-0.6,0.6)
    ax.set_ylim(-1.1,0.2)
    ax.set_aspect(1)
    ax.set_frame_on(False)
    ax.set_xticks([])
    ax.set_yticks([])
    fig = plt.gcf()
    fig.set_figheight(8)

Once that's defined, generating a single image is just a function call.

In [2]:
pic(30*np.pi/180)

You could then generate a sequence of frames using code like so:

In [3]:
thetas = np.linspace(0,20*np.pi/180,20)
thetas = np.concatenate([thetas,thetas[-2:0:-1],-thetas,-thetas[-2:0:-1]])

for i in range(len(thetas)):
    if(i<10):
        file_name = "pic00"+str(i)+".png"
    elif(10<=i<100):
        file_name = "pic0"+str(i)+".png"
    else:
        file_name = "pic"+str(i)+".png"
    pic(thetas[i]);
    fp = open(file_name,'w')
    plt.savefig(file_name)
    fp.close()
    plt.close()

Finally, you can assemble these PNGs into an animation using ImageMagick. You can execute this from the command line like so:

convert -delay 1 pic*.png -loop 0 anim.gif

Problem is, the result doesn't look quite right. To improve this, you need to solve the pendulum equation $$\theta'' = -c\sin(\theta)$$ with initial conditions $$\theta(0) = 0 \: \text{ and } \: \theta'(0) = s_0$$ and use the resulting values of $\theta$ to generate the animation. The constant $c$ depends on the length of the pendulum but we'll just assume that $c=1$. The constant $s_0$ is the initial angular speed of the pendulum. Note that $s_0$ is particular to you. My daughter would compute her $s_0$ via

seed('audrey')
s0 = 0.2+random()/4

In addition, you should compute the period of the pendulum using interp1d and brentq. You'll solve the equation over a large time interval (10 or 20 seconds) to generate a list of $\theta$ values. Then use scipy.interpolate.interp1d to generate an interpolation function for the thetas. To get the period, use scipy.optimize.brentq to find when $\theta(t)=0$. Finally, re-solve the equation IVP over the proper time interval.

So, I guess this looks a little tricky - certainly more free flowy than our previous projects. Note, though, that I showed you something similar the first day of class when we examined planetary motion. Also, I've opened a discussion on discourse that provides you the opportunity to fully discuss a similar problem.