Notes and Examples for Python
For Linux the packages numpy, scipy, Gnuplot, and matplotlib should be in usr/lib/Python2.x/site-packages. {I like Gnuplot better than matplotlib.}
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
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 *
I prefer Gnuplot.py for plotting.
Note, 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.
******
In Gnuplot.py ver 1.7 the mouse support in linux (Fedora 7) doesn't work. To fix it do these 3 things:
Comment out line 188 {self('set terminal %s' ...} in _Gnuplot.py.
******
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
Plotting with matplotlib
plot(x,y) #To make and show plots using matplotlib (matlab-like)
grid(b=1) # or grid(True).
show() # But I couldn't get programs to continue after show().
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.
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: 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')