iPython Guide

Computational methods from PHYS 250
–>Download the Anaconda distribution of Python at http://continuum.io/downloads

Week 1: Root Finding
import scipy.optimize as opt
(root, info) = opt.brentq(f,a,b,full_output=True)
iterations = info.iterations
*brentq may not always work, but is faster than bisect
**plot f to determine bounds a and b
***use axhline(y=0) to graphically find roots
****use axvline(x=root) to check your answer

Week 2: Interpolation
import scipy.interpolate as inter
sp = inter.InterpolatedUnivariateSpline(x,y)
abs_err = abs(true(val)-sp(val))
rel_err = abs(true(val)-sp(val))/true(val)
*extrapolation, or interpolation beyond the values contained in array x and/or array y, is dangerous!
**you can also use the lagrange function to interpolate

Week 3: Derivative
%loadpy http://www.phys.cwru.edu/courses/p250/examples/richardson_center.py
rich_array = richardson_center(f,val,step_size,nstep)
best_deriv_est = rich_array[-1,-1]
*making the step size smaller does not always reduce the error
**you can modify richardson_center to employ richardson forward or backward differencing techniques
***you can also use sp.derivatives to take derivatives of a spline
****you can use sp.derivative(1) to get the derivative as a function

Week 4: Integration
import scipy.integrate as integ
(res,err,info) = integ.quad(f,a,b,full_output=True)
neval = info[‘neval’]
*integ.newton_cotes returns the weights and error coefficient for Newton-Cotes integration
**integ.newton_cotes(1) refers to trapezoidal integration
***integ.newton_cotes(2) refers to Simpson’s integration
****you can also use integ.trapz(f(x),x) to get the integral of function f over x
*****integ.cumtrapz works like integ.trapz, but shows a table
******integ.simps works like integ.trapz, but with better accuracy
*******integ.romberg works like integ.quad, but setting show=True shows a table

Week 5 and 6: Solving Systems of Differential Equations
import scipy.integrate as integ
def f(y,t):
return dydt
res = integ.odeint(f,y0,t)
*f is oftentimes the force equation, which is directly related to acceleration
**y0 is the intial conditions of y
***t is an array of different times
****res is a 2D array of shape (len(t),len(y))

Week 8: Eigenvalue Problems
import scipy.linalg as la
(lam,B) = la.eig(A)
*only requires one input! 2D array A
**there are many different ways to determine if A is a singular matrix
***singular matrices make inverse calculations unreliable
****you can still compute eignvalues and eigenvectors accurately from a singular matrix

Week 7: Solving Systems of Equations with Linear Algebra
import scipy.linalg as la
res = la.solve(A,b)
*A is the array typically associated with the left-hand side of the systems of equations
**b is the array typically associated with the right-hand side of the systems of equations
***you can never get the wrong answer! check by dot(A,res)-b. you should get zeros
****you can determine whether you solved for a linear or non-linear matrix by la.lu(A)

Week 9: Singular Value Decomposition
import scipy.linalg as la
(U,S,VT) = la.svd(A)
Sinv = 1./S
xp = dot(VT.T,Sinv*dot(U.T,b))
z = VT[-1:].T
x = xp + sum(alpha*z,axis=1)
*x is the general solution to singular matrix A, so alpha is any random number
**you can also construct a psuedo inverse by using the where command or some equivalent short-hand method

Week 10: Chi Squared Fitting
import scipy.optimize as opt
import scipy.special as sf
def f(x,*paramters):
return something
(a1,C1) = opt.curve_fit(f,xdata,ydata,sigma=sig,p0=guess_array)
chisq1 = sum((ydata-f(x,*a1))**2/sig**2)
dof1 = len(ydata)-len(a1)
gof1 = sf.gammainc(dof1/2.,chisq1/2.)
*you can plot error bars by using the errorbar function
**sqrt(diag(C1)) is the error associated with the values stored in a1
***sigma=sig is an optional input if you have the error associated with ydata
****p0=guess_array is an optional input to troubleshoot bad multi-dimensional fits

Week 11: Fast Fourier Transform
Hn = fft.fft(hn)
freq = fft.fftfreq(N)
ind = arange(1,N/2.+1)
psd = abs(Hn[ind])**2+abs(Hn[-ind])**2)
plot(freq[ind],psd,label=’Power Spectral Density of hn’)
sigma = val
tt = arange(-N/2.,N/2.)
g = exp(-tt**2/2/sigma**2)/sqrt(2*pi)/sigma
hn_smooth = convolve(hn,g)
plot(t,hn_smooth)
*hn is a noisy signal constructed by N number of points
**you can find dominating frequencies with the power spectral density
***you can also convolute by hand
****convolution by computer results in zero-padding
*****you can mimic filters by using the where function or some other equivalent short-hand method

Week 12: Random Numbers
def f(x):
return something
def p(y):
return something
x = rand(N)
y = f(x)
bins = linspace(a,b,num)
hist(y,bins=bins,histtype=’step’,normed=True)
plot(bins,p(bins))
*you can never get the probability transformation wrong. the histogram should match the plot provided N is large enough

Example Problem 1: Solving x = 2^(-x)
def f(x):
return 2.**(-x)-x
root = opt.brentq(f,0,1)
x = linspace(0,1,1000)
plot(x,f(x),’r’,label=’f(x)’)
axhline(y=0)
axvline(x=root)
legend(loc=’best’)
#root = 0.641

figure_1

Example Problem 2: Estimate the Period of Satellite Orbit
import urllib
url = ‘ http://www.phys.cwru.edu/courses/p250/data/hw2.dat
(t,pos) = loadtxt(urllib.urlopen(url), unpack=True)
sp = inter.InterpolatedUnivariateSpline(t,pos)
root1 = opt.brentq(sp,0,0.5)
root2 = opt.brentq(sp,10,10.5)
x = linspace(-1,11,10000)
plot(t,pos,’r+’,label=’data’)
plot(x,sp(x),’r’,label=’interpolation’)
axhline(y=0)
axvline(x=root1)
axvline(x=root2)
xlim(-1,11)
legend(loc=’best’)
#period = abs(root2-root1) = 10.079

figure_2

Example Problem 3: Estimate the Period of Satellite Orbit
sp = inter.InterpolatedUnivariateSpline(t,pos)
sp_deriv = sp.derivative(1)
root3 = opt.brentq(sp_deriv,2,3)
root4 = opt.brentq(sp_deriv,7,8)
x = linspace(-1,11,10000)
plot(t,pos,’r+’,label=’data’)
plot(x,sp(x),’r’,label=’interpolation’)
plot(x,sp_deriv(x),’g’,label=’derivative’)
axhline(y=0)
axvline(x=root3)
axvline(x=root4)
xlim(-1,11)
legend(loc=’best’)
#period = 2*abs(root4-root3) = 10.396

figure_3

Example Problem 4: Estimate Radius of Satellite Orbit
half_circumference = integ.quad(sp_deriv,root3,root4)
circumference = 2*half_circumference[0]
#r = circumference/(2*pi) = 3.817

figure_4

Example Problem 5: Projectile Motion with Linear Air Resistance bv
b = 0.002
m = 0.11
g = 9.8
def projectile_force(y,t):
dydt = zeros_like(y)
dydt[0] = y[1]
dydt[1] = -g-b*y[1]/m
return dydt
t = linspace(0,3,1000)
y0 = [0,8]
res1 = integ.odeint(projectile_force, y0, t)
ind = where(res1[:,0] == max(res1[:,0]))
plot(t,res1[:,0],’r’,label=’Position’)
plot(t,res1[:,1],’g’,label=’Velocity’)
axhline(y=0)
title(‘Projectile Motion with Linear Air Resistance bv’)
xlabel(‘Time’)
xlim(0,t[ind])
ylim(-5,10)
legend(loc=’best’)

figure_5

Example Problem 6: Projectile Motion with Non-Linear Air Resistance cv^2
c = 0.002
def projectile_force(y,t):
dydt = zeros_like(y)
dydt[0] = y[1]
dydt[1] = -g-c*y[1]**2/m
return dydt
res2 = integ.odeint(projectile_force, y0, t)

figure_6

Example Problem 7: Solve for Eigenvalues and Eigenvectors of Matrix A
A = ([[4,-7,3],[1,3,-3],[3,-29,21]])
(lam,B) = la.eig(A)
#lam = [8.88e-16, 2.78e+00, 2.52e+01]
#since one of our eigenvalues are close to 0, we know A is a singular matrix
#can we trust our lamda computation?
#diag(dot(dot(inv(B),A),B)) = [-1.22e-15, 2.78e+00, 2.52e+01]
#we can reliably solve for the eigenvalues of matrix A

Example Problem 8: Solve the Systems of Equations Defined by Matrix A and b
b = ([-1,-2,8])
res = la.solve(A,b)
#res = [-0.579, 0.0263, 0.5]
#can we trust our solution?
#dot(A,res)-b = [0,0,0]
#this is a solution, but is this solution unique?
ludecomp = la.lu(A)
Umatrix = ludecomp[2]
#Umatrix[-1,-1] = -4.44e-16
#this is close to 0, so our solution is not unique. we have solved for a linearly-related matrix

Example Problem 9: Write the Complete Solution to the Systems of Equations Defined by Matrix A and b
(U, S, VT) = la.svd(A)
Sinv = 1./S
Sinv[2] = 0
xp = dot(VT.T,Sinv*dot(U.T,b))
z = VT[-1:].T
#x = xp + sum(alpha*z,axis=1)
#check to see if this works
alpha = rand(1)
x = xp + sum(alpha*z,axis=1)
#dot(A,x)-b = [-1.33e-15, -4.44e-16, 1.78e-15]

Example Problem 10: Fitting a Curve to a Gaussian Distribution
(y,b) = histogram(randn(10000),bins=50)
b_center = (b[:-1] + b[1:])/2
def f(x,A,u,s):
return A*exp(-0.5*((x-u)/s)**2)
(a1,C1) = opt.curve_fit(f,b_center,y)
#chisq = sum((y-f(b_center,*a1))**2/y) = 48.197
#dof = len(y)-len(a1) = 47
#gof = sf.gammaincc(dof/2.,chisq/2.) = 0.424
plot(b_center,y,’r’,label=’Gaussian distribution’)
plot(b_center,f(b_center,*a1),’g’,label=’Best fit’)
title(‘Best Fit of Random Gaussian Distribution’)
xlim(-3.5,3.5)
legend(loc=’best’)

figure_10

Example Problem 11: Smoothing a Noisy Signal
N = 10000
def f(x):
return sin(pi*x)
x = linspace(0,10,N)
hn = f(x)+randn(N)
sigma = 100.
tt = arange(-N/2,N/2)
g = exp(-tt**2/2/sigma**2)/sqrt(2*pi)/sigma
hn_smooth_zero_padded = convolve(hn,g)
hn_smooth = zeros(10000)
hn_smooth[:] = hn_smooth_zero_padded[5000:15000]
plot(x,hn,’r’,label=’Noisy Signal’)
plot(x,hn_smooth,’g’,label=’Smooth Signal’)
title(‘Smoothing of a Noisy Signal’)
xlabel(‘Time’)
ylabel(‘Amplitude’)
legend(loc=’best’)

figure_11

Example Problem 12: Monte Carlo Integration of 1/(sqrt(x)*(e^x+1))
def f(x):
return 1./(sqrt(x)*(e**x+1))
true = integ.quad(f,0,2)
def p(y):
return k/sqrt(y)
k = 1./(2*sqrt(2))
x = rand(1000000)
y = 2*x**2
ys = sort(y)
res = trapz(f(ys),ys)
rel_err = abs(true[0]-res)/true[0]
#rel_err = 8.058e-07
#res = 0.997
bins = linspace(0,2,1000)
hist(y,bins=bins,histtype=’step’,normed=True)
plot(bins,p(bins),’r’)
title(‘Transformation Confirmation’)
ylim(0,4)

figure_12

Advertisements

Feel free to comment

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s