Notes and Examples for Python
For Linux the packages numpy, scipy, Gnuplot, and matplotlib should be in usr/lib/Python2.x/site-packages.
Arrays, Numerical Integration, Numerical Solution of ODEs, Curve Fitting, Reading and Writing Array files, Finding zeros of functions, Graphing with Gnuplot, Fast Fourier Transform, Waveforms: Square, Sawtooth, Time Delay, Noise, Create Postscript Graph, Simple Plots with matplotlib, Interactive Plots with matplotlib, Plotting with log or linear axes
Get python from www.python.org, scipy (and numpy) from www.scipy.org, matplotlib from http://matplotlib.sourceforge.net/, gnuplot-py from http://sourceforge.net/project/showfiles.php?group_id=17434
Good references: http://www.scipy.org/Cookbook, for example, to learn plotting: http://www.scipy.org/Cookbook/Matplotlib, for gnuplot, see demo.py and test.py in the gnuplot-py directory
To learn about packages, subpackages, and routines, run python from a terminal, then follow these example commands to learn about the scipy package:
To make programs executable for Linux 1st line of program is:
#! /usr/bin/env python {and change permissions to include execute}
Programs should load packages as follows:
from __future__ import division
from scipy import *
from numpy import * # This not necessary since scipy loads numpy
from numpy.random import * # If need random number generator
import Gnuplot # If use gnuplot for graphs
from matplotlib import * # If use matplotlib for graphs
from pylab import *
Note, if using Gnuplot.py, then gnuplot must also be installed. http://www.gnuplot.info
Problems installing Gnuplot.py on Fedora
8 (I don't remember if I had problems with 7). Doing "python setup.py
install" did not work. So just copy the unzipped gnuplot-py-1.7 directory
to "site-packages" and rename it Gnuplot.
To enter values interactively use raw_input(text). (Doesn't work on emacs!! Use Eric)
x = float(raw_input('Enter x value: ')) #converts text input into real.
x = arange(0,10) # To make integer array vector of x values.
x = linspace(0,10,101) # To make more general array.
y = x**2
x=zeros(100,float) # creates x array filled with 0
Example: To include Time Delay
import time
time.sleep(2.0) #Causes 2.0 second delay.
Example: Numerical integration
from scipy import *
from scipy.integrate import quad
"""
Example of integration of a function using quad
Do the exponential function.
"""
result = quad(lambda x: exp(-x),0,5) #Uses anonymous function, lambda
disp(['Numerical result: ',result])
analytic = 1-exp(-5)
disp(['Analytic result: ',analytic])
Example: Numerical solution of ODEs
from scipy import *
from scipy.integrate import odeint
from matplotlib import *
from pylab import *
"""
Example of solving system of differential equations.
Do the damped oscillator, b is damping factor
"""
def damped_osc(u,t,b): #defines the system of odes
x=u[0]
v=u[1]
return(v,-x-b*v) #the derivatives of u
t = arange(0,20,0.1)
u0 = array([1,0])
b=0.4
u=odeint(damped_osc,u0,t,args=(b,)) #b is in tuple, needs comma
figure(1) #Using matplotlib
plot(t,u[:,0],t,u[:,1])
xlabel('Time')
title('Damped Oscillator')
figure(2)
plot(u[:,0],u[:,1])
title('Phase-space')
xlabel('Position')
ylabel('Velocity')
show()
from scipy import *
from matplotlib import *
from pylab import *
from scipy.optimize import leastsq
"""
Example of curve fitting for
a1*exp(-k1*t) + a2*exp(-k2*t)
"""
def dbexpl(t,p):
return(p[0]*exp(-p[1]*t) + p[2]*exp(-p[3]*t))
a1,a2 = 1.0, 1.0
k1,k2 = 0.05, 0.2
t=arange(0,100,0.1)
data = dbexpl(t,[a1,k1,a2,k2]) + 0.02*randn(len(t))
def residuals(p,data,t):
err = data - dbexpl(t,p)
return err
p0 = [0.5,1,0.5,1] # initial guesses
guessfit = dbexpl(t,p0)
pbest = leastsq(residuals,p0,args=(data,t),full_output=1)
bestparams = pbest[0]
cov_x = pbest[1]
print 'best fit parameters ',bestparams
print cov_x
datafit = dbexpl(t,bestparams)
plot(t,data,'x',t,datafit,'r',t,guessfit)
xlabel('Time')
title('Curve-fitting example')
grid(True)
show()
Example: File IO, Reading and writing arrays
from scipy import *
from scipy.io import write_array,read_array
"""
Example of file io, reading and writing array files
"""
x=arange(0,10)
y=exp(-0.2*x)*cos(2*pi*x/20.)
data=column_stack((x,y))
print data
write_array("test_file_io.txt",data,separator=' ',linesep='\n',precision=8)
A = read_array("test_file_io.txt")
print A
print "done"
Example: Simple Plots with matplotlib
from scipy import *
from pylab import *
"""
Shows how to put a couple of plots on a graph using matplotlib
"""
x = arange(0, 10, 0.5) #Define x-range: 0 to 10 with stepsize of 0.5
y = x**2 # Function x squared
z = 5*x # Function 5x
#plot using blue circles with line, and red x's no line
plot(x,y,'bo-',label='$x^2$',markersize=8.0,linewidth=2.0)
plot(x,z,'rx',label='$5x$',markersize=10.0,linewidth=2.0)
legend()
xlabel('x axis')
ylabel('y axis')
title('Simple Graph')
grid(True)
show()
Example: Interactive Plots with matplotlib
#!/usr/bin/env python
from scipy import *
import matplotlib
matplotlib.use('TkAgg') #Try this to work best in linux
from pylab import *
print """Example of interactive plotting with matplotlib's pyplot.
Won't work correctly on Windows using IDLE. Run from Command Prompt"""
x = linspace(0,1) #Defines array of x values from 0 to 1
f = 2
ion() #Turns interactive mode on.
ans=''
while ans != 'q':
y = sin(2*pi*f*x)
plot(x,y,'o')
draw() #Need this to work on linux?
ans=raw_input('f = '+repr(f)+'. Enter new f, or q: ')
if ans != 'q':
f = float(ans)
clf() #Clears the figure
ioff()
close() #Closes the figure
#Use following 2 lines instead of close() if want figure to remain.
#print 'Remember to close figure window.'
#show()
Example: Switching between log and linear axes on plots
#! /usr/bin/env python
from __future__ import division
from scipy import *
import matplotlib
matplotlib.use('TkAgg')
from pylab import *
print "Example of switching between log and linear axes"
RC = 1e-3
f=logspace(1,5) #good for logarithmic axes
ans=' '
ion()
while ans != 'q':
y=1/sqrt(1+(2*pi*f*RC)**2)
plot(f,y,'b+-')
grid(True)
draw()
ans=raw_input('Enter log, lin, RC, q: ')
if ans == 'log':
semilogx()
if ans == 'lin':
delaxes()
if ans =='RC':
RC=float(raw_input('Old RC = '+repr(RC)+'. Enter new RC: '))
clf()
close()
ioff()
Example: Finding zeros of nonlinear functions using fsolve
from scipy import *
from scipy.optimize import fsolve
from matplotlib import pylab
from pylab import *
"""
Use fsolve to find intersections of exp(x)-1 with cos(x)
Note that need different guesses to find different intersections
"""
def f(x):
return (exp(x)-1)-cos(x)
result = fsolve(f,-1.5)
print result
x = linspace(-5,1,101)
plot(x,exp(x)-1,x,cos(x))
grid(b=1)
show()
Example: Plotting using Gnuplot
#!/usr/bin/env python
from numpy import *
import Gnuplot
"""
Demonstrates how to use gnuplot for graphs. Plots 2 graphs with different styles.
"""
x = linspace(0,10,31)
y = x**2
y2 = 10*sin(pi*x)
g = Gnuplot.Gnuplot() #!! Won't be able to use 'with' in python 2.6?
d = Gnuplot.Data(x,y,title='squared', with='lp lt 1 pt 6')
d2=Gnuplot.Data(x,y2,title='sine',with='p pt 7 ps 4')
g('set grid')
g.xlabel('x axis')
g.ylabel('y axis')
g.plot(d,d2)
ans = raw_input('Enter f to create .png file, or Enter to quit ')
if ans == 'f':
g.hardcopy('filename.png',terminal = 'png')
#g.hardcopy('filename.ps',terminal='postscript',enhanced=1,color=1)
g.reset()
from __future__ import division
from scipy import *
import Gnuplot
"""
Test the fft routine. Add signals, and multiply signals.
"""
t=linspace(0,1,513)
# Use 2^N + 1
dt = (t[-1]-t[0])/(len(t)-1) # Maximum frequency is 1/dt ?
fmax = 1/(2*dt)
f1 = 80
f2 = 90
sig = 1 + sin(2*pi*f1*t) + 1 + sin(2*pi*f2*t) # sum of signals
#sig=(1+sin(2*pi*f1*t))*(1+sin(2*pi*f2*t)) # product of signals
npts=512 #Use some power of 2
ft = fft(sig,n=npts)
mgft=abs(ft)
df = fmax/float(npts/2)
f=linspace(0,fmax,npts/2+1)
print 'fmax = ',fmax,' df = ',df,' ','\n f = ',f[0:5]
g=Gnuplot.Gnuplot()
d=Gnuplot.Data(f,mgft[0:npts/2+1],with='l')
g.xlabel('frequency')
g.ylabel('fft magnitude')
g.plot(d)
raw_input('Press Enter ')
g.reset()
Example: Make a square or sawtooth wave using scipy.signal.waveforms
#! /usr/bin/env python
from scipy import *
from scipy.io import write_array
from scipy.signal.waveforms import square, sawtooth
import Gnuplot
# Make a square wave or sawtooth. Plot it, save image, save data file
t = linspace(0,0.01,501)
f = 1000
A = 3
#y = A*square(2*pi*f*t,duty=0.5)
y = A*sawtooth(2*pi*f*t,width=0.5)
g=Gnuplot.Gnuplot()
d =Gnuplot.Data(t,y,with='l')
g('set yrange[-3.5:3.5]')
g.plot(d)
ans = raw_input('Enter i to make image file, +\n'+
'd to make data file, \n or press enter ')
if ans == 'i':
filename = raw_input('Enter filename, include .png ')
g.hardcopy(filename,terminal='png')
if ans == 'd':
filename = raw_input('Enter filename for data file ')
data=column_stack((x,y))
write_array(filename,data,separator=' ',precision=8)
g.reset()
#! /usr/bin/env python
from numpy import *
from numpy.random import normal
import Gnuplot as G
"""
Make signal with gaussian noise
"""
npts=200
x=linspace(0,10,npts)
theory=5*sin(2*pi*x/2.0)
noise=normal(0,1,npts) #mean, std dev, num pts
sig=theory+noise
g=G.Gnuplot(debug=1)
pre=G.Data(x,theory,with='l')
dat=G.Data(x,sig,with='p')
g.plot(dat,pre)
raw_input('press enter')