# Python part 2, solving ordinary differential equations

## Numerical solution of differential equations.

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()

## Analysis of the dampled pendulum

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.