workshop_3.py

import numpy as np
import pylab as lab

# don't show all decimal places when printing
np.set_printoptions(precision=4)
 
# define a, b and c
a = np.ones([4,4])
print "a:\n", a
 
b = a * 2
print "\nb:\n", b
 
c = np.arange(16)
c = c.reshape(4, 4)
print "\nc:\n", c
print
 
# Check that the traces of a,b and c are 4, 8 and 30 respectively
print "The traces of a, b and c are 4, 8 and 30 respectively:",
print [4,8,30] == [np.trace(m) for m in (a,b,c)]
print
 
# Check that the determinants of a,b and c are all 0
import numpy.linalg as npla
print "The determinants of a,b and c are all zero:",
print all([npla.det(m)==0 for m in (a,b,c)])
 
cs = np.sqrt(c)
print "\nsqrt(c):"
print cs
print "\nDeterminant of sqrt(c): ",
print npla.det(cs)
 
print "\nInverse of sqrt(c):"
print npla.inv(cs)
 
print "\nEigenvalues of a and b:"
a_eval, a_evec = npla.eig(a)
b_eval, b_evec = npla.eig(b)
print "a: ", np.round(a_eval, 5)
print "b: ", np.round(b_eval, 5)

print "\nThe eigenvalues of b are 2x eigenvectors of a"
 
print "\nEigenvectors (to 4 d.p.):"
print "a: "
for v in a_evec.transpose():
    print np.round(v, 5)

print "b: "
for v in b_evec.transpose():
    print np.round(v, 5)

print "\nThe eigenvectors of a and b are the same"

print "\nCheck a.eigenvector = eigenvalue*a:"
for val, vec in zip(a_eval, a_evec.transpose()):
	print np.round(np.dot(a, vec),2), "=", np.round(val*vec,2)

print "\nnp.eye(5):"
print np.eye(5)
print "np.eye creates an identity matrix of size NxN"
print
 
# 3B - slicing arrays
print "Part 3B"
q = np.loadtxt('numpy_test.txt')
print "2) q is a %d row by %d col matrix" % np.shape(q)
print "3) q[3,2] is %r" % q[3,2]
print "4) np.shape(q[2:6,1:4]) as a %d row by %d col matrix" % np.shape(q[2:6,1:4])
print "5) the largest element of q[2:3, :3] is %r" % np.max(q[2:3, :3])
print "6) every other element of the first column of q"
print "Even elements:"
# print [x for i, x in enumerate(q[:,0]) if i%2==0]
print q[::2,0]
print "Odd elements:"
# print [x for i, x in enumerate(q[:,0]) if i%2==1]
print q[1::2,0]
print "7) bottom right 3x3 matrix:"
print q[-3:,-3:]
 
# 3C Lagrange Interpolation
def l_i(x, xs, i):
    """The ith Lagrange basis polynomial, Li(x)."""
    xi = xs[i]
    num, den = 1.0, 1.0
    for j, xj in enumerate(xs):
        if not j == i:
            num *= (x - xj)
            den *= (xi - xj)
    return num/den
 
def L(x, xs, ys):
    """Returns the lagrange interpolated value of x, given points Xs and Ys."""
    ps = zip(xs, ys)
    return sum(y*l_i(x, xs, i) for (i, (_,y)) in enumerate(ps))
 
def plot_lagrange(xs, ys):
    lxs = np.linspace(min(xs), max(xs), 100)
    lys = [L(x, xs, ys) for x in lxs]
    lab.plot(lxs, lys, label='lagrange')
    lab.plot(xs, ys, label='points')
    lab.legend()
    lab.show()
 
lag = np.loadtxt('lagrange_data.txt')
lab.title("Plot of lagrange_data.txt interpolated using lagrange polynomial")
plot_lagrange(lag[:,0],lag[:,1])
 
 
# 3D Runge's Problem
def f(x):
    return (1.0 / (1 + 25*x**2))
 
def runge_xs(n):
    return np.array([(2.0*i)/n - 1 for i in range(n+1)])
 
x_runge = runge_xs(10)
y_runge = f(x_runge)
lab.title('Runge\'s Problem: At n=10 there are oscillations at the ends of the interval')
plot_lagrange(x_runge, y_runge)
 
 
lab.title( "Runge's Problem: At larger values of n, oscillations grow")
lab.plot(runge_xs(100), f(runge_xs(100)), label='f(x)')
for n in range(5, 15, 3):
    xs = runge_xs(n)
    lys = [L(x, xs, f(xs)) for x in np.linspace(-1,1,100)]
    lab.plot(np.linspace(-1,1,100), lys, label='lagrange n=%d' % n)
    lab.legend()

lab.show()
 
# 3E Chebyshev nodes
def cheby_xs(n):
    return np.cos([((2*i-1)/(2.0*n))*np.pi for i in range(1, n+1)])
    
# Cheby spaced points do not show such large oscillations
# and fit f(x) with lower error given a larger n
lab.title('Runge\'s Problem: Chebyshev distributed points do not show such large oscillations (n=21)')
plot_lagrange(cheby_xs(21), f(cheby_xs(21)))