# Python III: Plotting and Differential Equations

## Plotting

### Basics

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

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.

• Lines can be drawn as '-', '--', ':', or '-.'
• Points can be represented by 'o', 'x', '+', '.', or 's' for square.
• Basic colors are 'b', 'r', 'g', 'k', 'c', 'm', 'y'.

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

clf()
plot(xdata, ydata, 'mo')
show()

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.

clf()
plot(xdata, ydata, 'mo', markersize=10)
xlim(-1,10)
ylim(-1,10)
show()

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
clf()
plot(xdata, ydata, 'mo', x, y, 'k--', markersize=10, linewidth=2)
show()

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=',')
plot(data[:6,1],data[:6,2],'o-')
plot(data[6:,1],data[6:,2],'x-')
legend(['15 feet per sec', '21.8 feet per sec'])
ylim(0,0.8)
xlim(0,25)

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

figure(2)
imshow(randn(20,30), interpolation='none')
colorbar()

### 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.

• Vector formats. These formats can be explicitly resized without any loss of quality in fonts or images. Unfortunately, this extra sophistication has lead them to have limited browser support and less widely available editing tools (Use Inkscape or Adobe illustrator, for example)

• .pdf (P)ortable (D)ocument (F)ormat graphics are vector graphics that can be easily included in LaTex documents compiled with PDF. Because they are vector-formated, they are resizable without any quality loss. But they don't work in web-pages.

• .svg (S)calable (V)ector (G)raphics are a vector format that is human-read-writeable and can be natively included in web pages because it is implemented as XML code browsers can read. Again, since it is a vector format, resizing does not effect quality. Not all features of SVG images are supported by all browsers.

• .eps (E)ncapsulated (P)ost(s)cript is a varient of the postscript printer language that draws vector graphics, can be edited by hand, and can be resized. EPS can be included in documents compiled with latex. But EPS is not well supported by windows, doesn't work in browsers, and is not widely used any more.

• Bitmat formats. These formats represent an image as a rectangular lattice of pixels with some set of numbers for each pixel describing it's color and transparency. These formats are simple and easy to render on a monitor, and so they gained traction quickly on the web and subsequently with digital photography. All are viewable in all browsers. However, they typically handle text badly under resizing, so they are not well-suited to diagrams.

• .jpg This common format uses lossy compression that reduces quality but makes file sizes very small, and hence cheap to store and send over internet.

• .gif GIF's are an old lossless bitmap format that used to have licensing issues, but is now free and has resurged in popularity because it can be used for simple animations also. Browser-native. Like all bitmap formats, resizing often leads to artifacts and quality loss.

• .tif TIF is a high-quality bitmap format that stores pixels explicitly and useful for images-as-data, like telescope images and OCR work. Because it of it's explicit representation, tif's are typically large files. Browser native.

• .png A newer, lossless compressed format but file sizes can be bigger than JPEG sizes.

• .mpg MPEG is a more powerful movie format for storing animations.

## 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
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*exp(r*timeEnd)
XFinalError = abs( XFinalCalc - XFinalExct )
errorStr = "Final Error = %1.5g"%(XFinalError)

figure(1)
options = {'fontsize': 18}

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

### 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$$.

• 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.

• The crazy double-pendulum 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.

## 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

### References

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

• $$x$$ is the angle of the governor,
• $$y$$ is the rate of change of the angle of the governor,
• $$z$$ is the rotation rate of the main drive.

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

• $$a$$ is the maximum drive force
• $$b$$ is the load force
• $$f$$ is the linear friction on the governor

#### Analysis

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.