Scientific Computing

Scientific computing involves the following:

Arrays, One dimensioned, Two dimensioned, Multi dimensioned
Numerical methods
Data visualization, plotting graphs in 2D and 3D

NumPy does not have array data type in the base language, instead it has List. Lists are versatile but inefficient for numerical computations. Required packages are:

NumPy - general n-dimensioned array implementation
SciPy - efficient numerical methods
Matplotlib - 2D plotting functions similar to graphing functions in MATLAB

Creating arrays

Arrays can be created from lists or other container types, such as tuples.

Function to create lists

range(n) - Create a list starting from zero, up to but not including n, incremeting at 1
range(n1, n2) - Create a list starting from n1, up to but not including n2, incrementing at 1
range(n1, n2, delta) - Create a list starting from n1, up to but not including n2, incrementing at delta, where delta is an integer which is positive if n1 < n2 and negative if n2 < n1.



In [1]:

    
a = range(10)
print a, type(a)









    



[0, 1, 2, 3, 4, 5, 6, 7, 8, 9] <type 'list'>



In [2]:

    
b = range(1, 10)
print b, type(b)









    



[1, 2, 3, 4, 5, 6, 7, 8, 9] <type 'list'>



In [3]:

    
c = range(2, 10, 2)
print c, type(c)









    



[2, 4, 6, 8] <type 'list'>



In [9]:

    
d = range(10, 0)
print d, type(d)









    



[] <type 'list'>



In [10]:

    
d = range(10, 0, -2)
print d, type(d)









    



[10, 8, 6, 4, 2] <type 'list'>

Creating Arrays

Arrays with equally spaced points

arange(n) - Similar to range() but returns an array instead of a list and delta can be a real number



In [4]:

    
from numpy import arange

a = arange(5)
print a, type(a)









    



[0 1 2 3 4] <type 'numpy.ndarray'>



In [5]:

    
b = arange(1, 5)
print b, type(b)









    



[1 2 3 4] <type 'numpy.ndarray'>



In [6]:

    
c = arange(1, 5, 0.5)
print c, type(c)









    



[ 1.   1.5  2.   2.5  3.   3.5  4.   4.5] <type 'numpy.ndarray'>



In [7]:

    
d = arange(10, 5)
print d, type(d)









    



[] <type 'numpy.ndarray'>



In [8]:

    
d = arange(5, 0, -0.5)
print d, type(d)









    



[ 5.   4.5  4.   3.5  3.   2.5  2.   1.5  1.   0.5] <type 'numpy.ndarray'>

Linearly spaced points

linspace(x1, x2, np) generates np data points starting with x1 up to and including x2. Note, np is the number of data points and therefore the number of equal intervals is n-1.



In [9]:

    
from numpy import linspace

a = linspace(0, 5, 11)
print a, type(a), a.ndim









    



[ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5  5. ] <type 'numpy.ndarray'> 1



In [10]:

    
b = linspace(5, 0, 11)
print b, type(b), b.ndim









    



[ 5.   4.5  4.   3.5  3.   2.5  2.   1.5  1.   0.5  0. ] <type 'numpy.ndarray'> 1



In [11]:

    
from numpy import pi

x = linspace(0, 2*pi, 501)
print type(x), x.ndim, x.shape, x.size









    



<type 'numpy.ndarray'> 1 (501,) 501

Explicit creation of arrays

NumPy provides the function array() to create arrays from lists. Elements of the list are used to populate the elements of the array.



In [12]:

    
from numpy import array

a = array([1, 2, 3, 4, 5])
print a, type(a), a.ndim, a.shape, a.size









    



[1 2 3 4 5] <type 'numpy.ndarray'> 1 (5,) 5

To create a two dimensioned array, use a two dimensioned list



In [13]:

    
b = array([ [1, 2, 3], [4, 5, 6] ])
print b, type(b), b.ndim, b.shape, b.size









    



[[1 2 3]
 [4 5 6]] <type 'numpy.ndarray'> 2 (2, 3) 6

Quick arrays

zeros() - array with all elements zeros
ones() - array with all elements ones
eye() - Identity matrix
diag() - Matrix with elements in given list on main diagonal or parallel to main diagonal

When creating arrays this way, keep in mind the following points:

Size of the array must be specified as a tuple
Data type of elements of the array must be specified



In [14]:

    
import numpy as np

a = np.zeros((3,5), dtype=float)
print a, type(a), a.ndim, a.shape, a.size









    



[[ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]] <type 'numpy.ndarray'> 2 (3, 5) 15



In [15]:

    
b = np.ones((5,3), dtype=int)
print b, type(b), b.ndim, b.shape, b.size









    



[[1 1 1]
 [1 1 1]
 [1 1 1]
 [1 1 1]
 [1 1 1]] <type 'numpy.ndarray'> 2 (5, 3) 15



In [16]:

    
c = np.eye(5, dtype=float)
print c, type(c), c.ndim, c.shape, c.size









    



[[ 1.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.]
 [ 0.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  1.]] <type 'numpy.ndarray'> 2 (5, 5) 25



In [17]:

    
d = np.diag([10, 20, 30])
print d, type(d), d.ndim, d.shape, d.size, d.dtype









    



[[10  0  0]
 [ 0 20  0]
 [ 0  0 30]] <type 'numpy.ndarray'> 2 (3, 3) 9 int64



In [18]:

    
x = np.diag([1.0, 2.0, 3.0])
print x, type(x), x.ndim, x.shape, x.size, x.dtype









    



[[ 1.  0.  0.]
 [ 0.  2.  0.]
 [ 0.  0.  3.]] <type 'numpy.ndarray'> 2 (3, 3) 9 float64



In [19]:

    
y = np.diag([10, 20, 30], 1)
print y









    



[[ 0 10  0  0]
 [ 0  0 20  0]
 [ 0  0  0 30]
 [ 0  0  0  0]]



In [20]:

    
z = np.diag([10, 20, 30], -1)
print z









    



[[ 0  0  0  0]
 [10  0  0  0]
 [ 0 20  0  0]
 [ 0  0 30  0]]



In [21]:

    
a = np.random.rand(3,5)
print a









    



[[ 0.13861321  0.78001907  0.72302429  0.50430888  0.17128537]
 [ 0.62910545  0.58665316  0.62767661  0.65094655  0.59628018]
 [ 0.00581169  0.05066921  0.95133704  0.66546244  0.04972209]]



In [22]:

    
b = np.random.rand(5, 3) * 10 + 1
print b









    



[[  3.28465335   5.76424024   5.76067066]
 [  3.51566546  10.98832049   3.23687704]
 [  2.5757205    9.02249089   1.91224102]
 [ 10.94410301   9.80179106   1.92447814]
 [  5.13211609   1.06601605   3.68179162]]

Common Operations on Arrays



In [23]:

    
a = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9]])
b = np.array([ [2, 4, 6], [8, 10, 12], [14, 16, 18] ])
c = a + b
print b









    



[[ 2  4  6]
 [ 8 10 12]
 [14 16 18]]



In [24]:

    
d = a * b
print d









    



[[  2   8  18]
 [ 32  50  72]
 [ 98 128 162]]



In [25]:

    
x = 2 * a
print x









    



[[ 2  4  6]
 [ 8 10 12]
 [14 16 18]]



In [26]:

    
y = 10 * a - 2 * b
print y









    



[[ 6 12 18]
 [24 30 36]
 [42 48 54]]



In [27]:

    
from matplotlib import pyplot as plt
%matplotlib inline

x = np.linspace(0, 2*np.pi, 501)
y1 = np.sin(x)
y2 = np.cos(x)

plt.plot(x, y1, label='sin(x)')
plt.plot(x, y2, label='cos(x)')
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Harmonic function')
plt.legend(loc=3)
plt.show()

Linear Algebra

Solve linear algebraic equations
Find determinant, eigenvalues, eigenvectors
LU decomposition
Norm, rank



In [28]:

    
a = np.array([ [10, 3, -4], [2, 8, -1], [3, 5, -9] ], dtype=float)
b = np.array([25, -20, 30], dtype=float)
x = np.linalg.solve(a, b)
print x









    



[ 1.66080844 -3.50615114 -4.72759227]



In [29]:

    
print np.dot(a, x) - b









    



[  0.00000000e+00   3.55271368e-15   0.00000000e+00]



In [30]:

    
print np.linalg.det(a)



In [31]:

    
print np.linalg.norm(a)









    



17.5783958312



In [32]:

    
w, X = np.linalg.eig(a)
print w
print X









    



[ 10.65003706   6.53077508  -8.18081215]
[[-0.81556075 -0.52860832  0.20930797]
 [-0.51871719  0.83249528  0.03452387]
 [-0.25650172  0.16590617  0.97724013]]

References

NumPy - NumPy home page
SciPy - SciPy home page
Pandas - Python data analysis library home page
SciPy Lecture Notes - Free book
Quantitative Economics - Free book