The simplest autonomous ordinary differential equation of interest in science is the exponential decay equation

\[ \begin{gather*} dn/dt = r n \quad \text{ where $r < 0$.} \end{gather*} \]

Implement the following script

```
# Numerical integration of the exponential decay model
#
# dn/dt = r n, where r < 0
#
from numpy import linspace, array
import scipy.integrate
from scipy import exp
from matplotlib.pyplot import \
plot, figure, text, show, savefig, xlabel, ylabel
r = -0.05 # decay_rate per day
def vector_field(X, t):
n = X[0]
return array([ r * n ])
time_start = 0.
time_end = 10.
numsteps = 20
observation_times = linspace(time_start, time_end, numsteps)
X_initial = array([32.])
X = scipy.integrate.odeint(vector_field, X_initial, observation_times)
X_final_calc = X[-1, 0]
X_final_exct = X_initial[0]*exp(r*time_end)
X_final_error = abs( X_final_calc - X_final_exct )
error_str = "Final Error = %1.5g"%(X_final_error)
figure(1)
options = {'fontsize': 18}
plot(observation_times, X[:, 0],'x-')
xlabel('Time (days)',options)
ylabel('Mass (micrograms)',options)
text(2,20,error_str,options)
savefig('exponential.pdf')
show()
```

A pendulum on a rod, hanging from a pivot with weak friction should move according to the following system of equations, according to Newton's laws.

\[ \begin{align*} \frac{d \theta}{dt} &= \omega, \\ \frac{d \omega}{dt} &= -\sin(\theta) - a \omega. \end{align*} \]

The angle is measured as different from the rest position. In our case, let's assume the friction coefficient is \(a=0.05\).

If we start with the pendulum not moving but at initial angle \(\theta = 3 \pi / 4\), make a plot showing the angle as a function of time for 3 full periods of the oscillation. How long does it take? Include a second subplot under the first that shows the angular velocity \(\omega(t)\) using matplotlib's subplot command. Make sure to put appropriate axes labels in your plot.

Download and run this demo of matplotlib's streamplot function for drawing phase planes.

Modify the streamplot demo to generate a phase plane of our pendulum equation.

At lunch, Wednesday, Shuhao was talking about the crazy double-pendulum, whose behavior takes allot of work to understand. Download and run this python animation to see it in action. Then, change the script so it illustrates our single-pendulum.

Starting at the same angle (\(\theta = 3 \pi / 4\)), find the initial velocity that lead to pendulum to stop exactly at the upside-down position after 1 full rotation. You can do this by trial and error, but I suggest instead using one of the root-finding methods like "ridder" from the last lab. (WARNING: this is a harder problem, so you should skip over it in class). Make an animation that illustrates your solution, and save it as a movie.