Python III: Plotting and Differential Equations




Python also has a useful set of plotting tools. The most commonly used and basic tools for plotting are from the matplotlib.pyplot package. There is an excellent gallery of examples. For a simple place to start, let us plot your solution to the least-squares problem, allong with the data, using Canopy's editor window. (use the m and b you calculated from the least-squares problem in lab II)

from matplotlib.pyplot import *

xdata = [0, 9, 5, 6]
ydata = [0, 9, 7, 5]
plot(xdata, ydata)

This doesn't really give us what we want -- we'd like points rather than a line. The plot function uses the same stylies for drawing lines that Octave and Matlab do.

So, to draw points as magenta circles with no connecting line, first we clear the figure using the clf() command. Then we redraw it.

plot(xdata, ydata, 'mo')

Still, all four points are not clearly visible -- they are small and stuck in the corners. We can fix this a little by manually elarging the limits on the domain and range of the plot, and by increasing the size of the points.

plot(xdata, ydata, 'mo', markersize=10)

Now, that looks good. Let's draw in your regression line now to see how it fits the data. (remember to fill in your values for m and b.

m = ?
b = ?
from numpy import linspace
x = linspace(0,10,3)
y = m * x + b
plot(xdata, ydata, 'mo', x, y, 'k--', markersize=10, linewidth=2)

Now, add title, an x-label, and a y-label to your plot. Make sure the font sizes are large enough to be easily read, even if the figure is shrunk. Save your figure to a pdf image file with the savefig command.

One extra useful hint. Often the data you want to plot are in point-form like pts = [ (0,0), (9,9), (5,7), (6,5) ] rather than the x-vector/y-vector form used by plot. The fastest way to convert from point-form to vector-form is with the zip command.

pts = [ (0,0), (9,9), (5,7), (6,5) ]
xdata, ydata = zip( * pts)

We can also plot data we've stored in a file. For example, this early data collected by George Cayley on lift in relation to the angle of attach of a wing.

figure() # to create a NEW figure
data = loadtxt('CayleyLiftData.csv', delimiter=',')
legend(['15 feet per sec', '21.8 feet per sec'])

There are other styles of plots as well. For example, you can plot images with imshow.

imshow(randn(20,30), interpolation='none')

Types of images

There are allot of different image formats you can choose from for saving your plots. Each has it's own good and bad sides.

Ordinary differential equation methods

Differential equations are some of the most useful ways of scribing natural systems and lead to some very interesting insights. The next examples show how to solve and analyze Ordinary differential equation systems with python.

Example 1: Exponential decay

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

\[ \begin{gather*} dn/dt = r n \quad \text{ where } \quad 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 vectorField(X, t):
    n = X[0]
    return array([ r * n ])

timeStart = 0.
timeEnd = 10.
numsteps = 20

observationTimes = linspace(timeStart, timeEnd, numsteps)

Xinitial = array([32.])

X = scipy.integrate.odeint(vectorField, Xinitial, observationTimes)

XFinalCalc = X[-1, 0]
XFinalExct = Xinitial[0]*exp(r*timeEnd)
XFinalError = abs( XFinalCalc - XFinalExct )
errorStr = "Final Error = %1.5g"%(XFinalError)

options = {'fontsize': 18}

plot(observationTimes, X[:, 0], 'x-')
xlabel('Time (days)', options)
ylabel('Mass (micrograms)', options)
text(2, 20, errorStr, options)

Example 2: The damped 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. \[ \frac{d^2\theta}{dt^2} = -\sin \theta - a \frac{d\theta}{dt}. \] This second-order equation can be rewritten as a system of two first-order equations by using the angular velocity \(\omega = d\theta/dt\) as a second variable. \[ \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\).

Example 3: Autocatalytic reactions

Once upon a time, it was believed that chemical reactions would always got to equilibrium straight-away, but some visionary chemists eventually showed that, in fact, chemical reactions can oscillate. The simpliest mathematical model of this phenomena is given by the bi-molecular mass-action reaction equations \[ \begin{align*} \frac{dp}{dt} =& a - p + p^2 q \\ \frac{dq}{dt} =& b - p^2 q \end{align*} \] where \(a\) and \(b\) are positive constants.

  1. Determine the steady-states of this system

  2. Find values of \(a\) and \(b\) for which the system is unstable and oscillates. Show with a time-series plot.

  3. Find different values of \(a\) and \(b\) for which the system is stable and converges to equilibrium. Show with a time-series plot.

Example 4: Watt's governor

Diagram of Watt's centrifugal flying-ball governor

Diagram of Watt's centrifugal flying-ball governor


Watt's governor was the first automatic control system used at the beginning of the industrial revolution for steam engines. It opened up whole new fields of both engineering and dynamical system analysis. In particular, James Maxwell showed that these governors could become unstable. Here are there equations, which we'll study a little.

\[ \begin{align*} \frac{dx}{dt} =& y \\ \frac{dy}{dt} =& z^{2} \sin{\left (x \right )} \cos{\left (x \right )} - \sin{\left (x \right )} - f y \\ \frac{dz}{dt} =& a \left(\cos{\left (x \right )}-b\right) \\ \text{parameter values} \quad & a = 1, \quad b = 1/2 \\ \text{initial condition} \quad & x(0) = 1.0472, \quad y(0) = 0., \quad z(0) = 1.386 \end{align*} \] Here, the state \((x, y, z) \in [0,\pi] \times (-\infty, \infty) \times [0, \infty)\). where

The parameters \((a, b, f) \in [0,\infty) \times [0,a] \times [0,\infty)\). where


  1. I. A. Vyshnegradskii and J. C. Maxwell independently showed that if the friction \(f\) is large enough, then the governor settles on a stable control. Plot timeseries for x,y,z when \(f = 4/5\) over 100 seconds.

  2. But if f is small enough, then the governor "hunts" for the best level without ever settling! Plot timeseries for x,y,z when \(f = 1/2\), and describe how the behavior differs from the previous case.

  3. Enter the vector field system into a column vector in sympy. Use sympy to solve for the steady-states of the system. Then, numerically estimate the steady-states for \((a,b,f) = (1,1/2,4/5)\) and \((a,b,f) = (1,1/2,1/2)\).

  4. Create a new function called "jacobian" to calculate the Jacob's derivative matrix of a vector field (given by a sequence of equations es) in terms of some variables (give by a sequence xs)

    >>> jacobian = lambda es, xs : Matrix([[ e.diff(x) for x in xs] for e in es])

    Use your new function to calculate the Jacobian of vector field for Watt's governor.

  5. Using the results of your calculates so far, create a new python the determines the sign of the largest eigenvalue of the Jacobian when evaluated at the steady-state solution. Numerically estimate the steady-states eigenvalues at steady-state for parameter sets \((a,b,f) = (1,1/2,4/5)\) and \((a,b,f) = (1,1/2,1/2)\). In which case is there an exponentially-growing part to solutions near the steady-state? At what value of \(f\) does this eigenvalue switch from having positive real part instead of a negative real part?

Example 5: Behavior Change in an epidemic

A simple extension of classic epidemic theory is a case where resistance to disease is behavioral, and individuals who become sick tend to abandon their protective behaviors. This kind of phenomena can be described with the system \[ \begin{align} \frac{dS}{dt} &= \mu - \beta S I + f \gamma I + a R - v S - \mu S, \\ \frac{dI}{dt} &= \beta (S+\sigma R) I - \gamma I - \mu I, \\ \frac{dR}{dt} &= - \sigma \beta R I + (1-f) \gamma I - a R + v S - \mu R. \end{align} \] where \(S\) is the fraction of susceptibles, \(I\) is the fraction of infecteds, and \(R\) is the fraction of recovered people (\(1=S+I+R\)). Calculate all steady-states when \(f=9/10\), \(\gamma=61\), \(a=1/6\), \(\mu =1/78\), \(\beta = 570\), \(\nu = 1/3\), and determine their local stability.