Although coding with Python is very versatile and allows many advanced features that are useful when manipulating massive data (a common task in science), Python is still a multipurpose language, what implies that scientific routines and functions cannot (should not) be supported within its basic core. Nevertheless, there are many different scientific libraries that can extend the capabilities of Python to scientific implementations in a natural way. Two of the most used libraries are NumPy and SciPy, intended for manipulating mathematical objects more efficiently and for extending and including numerical methods, respectively. Another less used libraries like SymPy are intended for manipulating analytical expressions, i.e. a CAS (Computer Algebraic System).
Installation of these libraries is often an easy task. In most of the Linux distros you should find them in the official repositories.
See the official pages of the libraries for new versions, news and manuals.
NumPy:
SciPy
SymPy
http://www.sympy.org/en/index.html
Anaconda
https://store.continuum.io/cshop/anaconda/
Anaconda is an interesting proposal that integrates many standard scientific libraries with Python
There are many different scientific libraries for Python with many different uses, even for very specific tasks. However, as we are interested in general numerical methods, we will focus only on NumPy and Scipy.
NumPy is the fundamental package for scientific computing with Python. It contains among other things:
Besides its obvious scientific uses, NumPy can also be used as an efficient multi-dimensional container of generic data. Arbitrary data-types can be defined. This allows NumPy to seamlessly and speedily integrate with a wide variety of databases.
NumPy can be imported in several different ways:
#Importing all the functions of NumPy
from numpy import *
#Importing NumPy as numpy
import numpy
#Importing NumPy with the alias of np (recommended)
import numpy as np
If imported with an alias, NumPy functions and features are accessed by using the alias before
#Basic functions
print np.exp(10.5), np.log(52.3), np.log10(63.9), np.sqrt(10.0)
36315.5026742 3.95699637107 1.80550085816 3.16227766017
#Trigonometric functions
print np.sin(5.0), np.cos(9.6), np.arcsin(0.5), np.arctan(5)
-0.958924274663 -0.984687855794 0.523598775598 1.37340076695
print "The value of PI is %1.6f"%( np.pi )
The value of PI is 3.141593
Supported methods for lists
x.append x.count x.extend x.index x.insert x.pop x.remove x.reverse x.sort
Supported methods for NumPy arrays
x.T x.clip x.dot x.item x.prod x.setfield x.take
x.all x.compress x.dtype x.itemset x.ptp x.setflags x.tofile
x.any x.conj x.dump x.itemsize x.put x.shape x.tolist
x.argmax x.conjugate x.dumps x.max x.ravel x.size x.tostring
x.argmin x.copy x.fill x.mean x.real x.sort x.trace
x.argsort x.ctypes x.flags x.min x.repeat x.squeeze x.transpose
x.astype x.cumprod x.flat x.nbytes x.reshape x.std x.var
x.base x.cumsum x.flatten x.ndim x.resize x.strides x.view
x.byteswap x.data x.getfield x.newbyteorder x.round x.sum
x.choose x.diagonal x.imag x.nonzero x.searchsorted x.swapaxes
Common
#Lists and numpy arrays can both store any type of data
x1 = [1.2, 3.5, 1.9]
x2 = np.array([1.6, -2.6, 6.9])
print x1, x2
[1.2, 3.5, 1.9] [ 1.6 -2.6 6.9]
#For lists, new elements can be added using append method
x1 = [1.2, 3.5, 1.9]
x1.append(5.9)
#For arrays, new elements can be added using append function of NumPy
x2 = np.array([1.6, -2.6, 6.9])
x2 = np.append(x2,3)
print x1,x2
[1.2, 3.5, 1.9, 5.9] [ 1.6 -2.6 6.9 3. ]
#A list can be converted into a numpy array
x = [1.1,3.4,1.0]
x = np.array(x)
#And a numpy array into a list as well
x = list(x)
Lists
#Operator + for lists is overloaded for concatenating
x1 = [1,2,3]
x2 = [3,2,1]
print x1+x2
[1, 2, 3, 3, 2, 1]
#Operator + for numpy arrays is overloaded for adding
x1 = np.array([1,2,3])
x2 = np.array([3,2,1])
print x1+x2
[4 4 4]
#Lists do not support other operators
x1 = [1.2,4.8,6.9]
x2 = [2.6,2.8,1.1]
#Multiplication
print x1*x2
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-16-10bb95b98344> in <module>() ----> 1 print x1*x2 TypeError: can't multiply sequence by non-int of type 'list'
#Division
print x1/x2
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-19-d643f42defa1> in <module>() 1 #Division ----> 2 print x1/x2 TypeError: unsupported operand type(s) for /: 'list' and 'list'
#Subtraction
print x1-x2
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-17-06680cf9cdd9> in <module>() ----> 1 print x1-x2 TypeError: unsupported operand type(s) for -: 'list' and 'list'
#Power
print x1**x2
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-18-da701bb7a380> in <module>() ----> 1 print x1**x2 TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'list'
NumPy arrays
#Numpy arrays support any mathematical operation (element by element)
x1 = np.array([1.2,4.8,6.9])
x2 = np.array([2.6,2.8,1.1])
print "Adding", x1+x2
print "Multiplication", x1*x2
print "Division", x1/x2
print "Subtraction", x1-x2
print "Power", x1**x2
Adding [ 3.8 7.6 8. ] Multiplication [ 3.12 13.44 7.59] Division [ 0.46153846 1.71428571 6.27272727] Subtraction [-1.4 2. 5.8] Power [ 1.6064649 80.81192733 8.37016462]
Note: Matrices can be represented as NumPy arrays where each element is a row vector. Nevertheless, be careful when multiply arrays, the operator * is overloaded in such a way that single elements are multiplied one by one, what is different than multiplication of matrices.
A = np.array([[1,2],[3,4]])
B = np.array([[4,3],[2,1]])
print A
print B
print A*B
[[1 2] [3 4]] [[4 3] [2 1]] [[4 6] [6 4]]
NumPy arrays have a number of interesting features that facilitate complex tasks. Next it is shown some of the more useful ones.
#It is possible to access elements of a numpy array using booleans
x = np.array([1,2,3,4])
y = np.array([True, False, True, False])
x[y]
array([1, 3])
#Operators >, <, >=, <= and == are also overloaded for numpy arrays
x = np.array([0,5,8,0])
y = np.array([0,6,5,1])
print x>y
print x<y
print x==y
print x>4
[False False True False] [False True False True] [ True False False False] [False True True False]
#Combining these features, we can perform searches and comparisons far more efficient
x = np.array([1,4,2,6,8,4,3,0,9,1,3,6,7,])
#A new list with numbers greater than 4
print x[x>4]
[6 8 9 6 7]
#Native methods of numpy arrays allow to calculate basic quantities
x = np.array([1,4,2,6,8,4,3,0,9,1,3,6,7,])
#Maximum element
print "Maximum element", x.max()
#Minimum element
print "Minimum element", x.min()
#Sorted arguments of the array
print "Sorted arguments", x.argsort()
#Sorted array
print "Sorted array", x[x.argsort()]
#Mean value
print "Mean value", x.mean()
Maximum element 9 Minimum element 0 Sorted arguments [ 7 0 9 2 6 10 1 5 3 11 12 4 8] Sorted array [0 1 1 2 3 3 4 4 6 6 7 8 9] Mean value 4.15384615385
NumPy includes many general purpose functions that complement the capabilities of python. We are interested here specially in functions for creating ordered arrays, storing and loading data as well as histograms, tasks that will be continuously required for the activities of the course.
#Create an array of 1's with a given size (even 2D sizes)
x = np.ones(5)
print x
[ 1. 1. 1. 1. 1.]
#Create an array of zeros with a given size (even 2D sizes)
x = np.zeros( (2,5) )
print x
[[ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.]]
#Create an array with a given range and a number of intervals
x = np.linspace( -np.pi, np.pi, 10 )
print x
[-3.14159265 -2.44346095 -1.74532925 -1.04719755 -0.34906585 0.34906585 1.04719755 1.74532925 2.44346095 3.14159265]
#Create an array with a given range and a given step
x = np.arange( 1, 5, 0.2 )
print x
[ 1. 1.2 1.4 1.6 1.8 2. 2.2 2.4 2.6 2.8 3. 3.2 3.4 3.6 3.8 4. 4.2 4.4 4.6 4.8]
#Using the function savetxt, it is possible to store data from a numpy array
data = np.array([[3.2, 2.1],[3.1, 4.1]])
np.savetxt( "file.dat", data, fmt="%1.5e %1.5e" )
#In the same way, using the function loadtxt it is possible to load external data files
data = np.loadtxt("file.dat")
#Data is then a multidimensional array with the loaded data
print data
[[ 3.2 2.1] [ 3.1 4.1]]
#Histograms are also useful and NumPy provides functions to do so
x = [1,5,3,7,2,4,6,8,9,5,2,3,5,6,7,8,3,4,5]
np.histogram(x)
(array([1, 2, 3, 2, 0, 4, 2, 2, 2, 1]), array([ 1. , 1.8, 2.6, 3.4, 4.2, 5. , 5.8, 6.6, 7.4, 8.2, 9. ]))
Other useful functions will be covered when needed during the course.
SciPy is a collection of mathematical algorithms and convenience functions built on the Numpy extension of Python. It adds significant power to the interactive Python session by providing the user with high-level commands and classes for manipulating and visualizing data. With SciPy an interactive Python session becomes a data-processing and system-prototyping environment rivaling sytems such as MATLAB, IDL, Octave, R-Lab, and SciLab.
Some of the packages included with SciPy are:
Each of these packages must be imported separately.
#Importing integrate package
import scipy.integrate as integ
The integrate package then includes the next functions:
integ.Tester integ.fixed_quad integ.odepack integ.quadrature integ.test
integ.complex_ode integ.newton_cotes integ.quad integ.romb integ.tplquad
integ.cumtrapz integ.ode integ.quad_explain integ.romberg integ.trapz
integ.dblquad integ.odeint integ.quadpack integ.simps integ.vode
Almost each of the numerical methods that will be covered during the course can be found in SciPy. In next classes we will explore the offered options by SciPy according to the specific methods covered.