Introduction to Numerical Methods/Python Programming

Python Programming

Objectives:

• familiarize with the Python programming language and command-line interface
  • use correct control structures
  • write simple Python programs with correct syntax
• import and use libraries

Resources:

• Introduction to Python for Computational Science and Engineering (A beginner’s guide)
• Google's Python class for developers
• The Beginner's Python Tutorial
• The Official Python Tutorial
• The Official Python for Beginners
• Introduction to Python
• IPython

Python is a scripting language. You can use the interactive Python console to write and execute code. After you log into a Unix/Linux system you can type the type the word "python" (with the quotes) to start the console and hit ^D (ctrl+D) will exit the console. This kind of interactivity makes it ideal for exploratory discovery - finding a solution by exploring and experimentation.

REPL

The Python console executes a Read-Eval-Print Loop (REPL), which means it prompts you for a expression/command and once you hit enter it will read your expression, evaluate it, print the value of the expression, and prompt you for the next expression. This makes it very easy to test a Python statement. You can use it as a calculator:

>>>1+2
3
>>>2**3
8
>>> # this is a comment
...
>>> 2**0.5
1.4142135623730951

Import Libraries

Many useful functions have been implemented in libraries. To use the functions you need to import them so that their names are recognized. The dir() function prints a directory of objects (functions) available in a module:

>>>import math
>>>dir(math)
['__doc__', '__name__', '__package__', 'acos', 'acosh', 'asin', 'asinh', 'atan', 'atan2', 'atanh', 'ceil', 
'copysign', 'cos', 'cosh', 'degrees', 'e', 'erf', 'erfc', 'exp', 'expm1', 'fabs', 'factorial', 'floor', 'fmod', 
'frexp', 'fsum', 'gamma', 'hypot', 'isinf', 'isnan', 'ldexp', 'lgamma', 'log', 'log10', 'log1p', 'modf', 
'pi', 'pow', 'radians', 'sin', 'sinh', 'sqrt', 'tan', 'tanh', 'trunc']
>>> help(math.exp)
Help on built-in function exp in module math:
exp(...)
    exp(x)
    
    Return e raised to the power of x.
>>> # this is similar to 2 ** 0.5
>>> math.sqrt(2)
1.4142135623730951

The help() function can be used to find the meaning and the usage of an object from a library as shown in the previous code example. Note that you need to hit the letter Q to quit the help session to return to the Python console.

Data Types

Python is not a strongly typed language, which means you do not need to declare a specific data type for a variable. But data items associated with variables have types, which are implied - derived from the expression or the operation on the data items. Therefore the data type of a variable can change over time. The type() function will tell you the type of data stored in a variable.

>>> type(a)
<type 'int'>
>>> type(1)
<type 'int'>
>>> a = 1.0
>>> type(a)
<type 'float'>
>>> a = "one"
>>> type(a)
<type 'str'>
>>> a = False 
>>> print a False
>>> type(a)
<type ’bool’>
>>> a = 1 + 2j
>>> type(a)
<type 'complex'>
>>> a = [1, 2, 3]
>>> type(a)
<type 'list'>

There are functions that convert between data types, e.g. int(), float(), str(), and list().

>>> a = int("123")
>>> a
123
>>> type(a)
<type 'int'>
>>> a = float("123.4")
>>> a
123.4
>>> type(a)
<type 'float'>
>>> a = str(123.4)
>>> a
'123.4'
>>> type(a)
<type 'str'>
>>> a =  list("one two three")
>>> a
['o', 'n', 'e', ' ', 't', 'w', 'o', ' ', 't', 'h', 'r', 'e', 'e']
>>> type(a)
<type 'list'>

Data types get promoted automatically as necessary.

>>> import sys
>>> sys.maxint
9223372036854775807
>>> type(sys.maxint)
<type 'int'>
>>> type(sys.maxint+1)
<type 'long'>

In Python strings, lists and tuples are sequences. They can be indexed and sliced in the same way. A list in Python can contain data items of different type.

>>> range(10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> for i in range(10):
...     print i ** 2
... 
0
1
4
9
16
25
36
49
64
81

Source Files

In addition to write and running Python code directly in the console you can save your commands/statements in a file and run the whole file as a script using the Python interpreter. For example, in a shell command-line you can type the following to run your program/script called hello.py:

python hello.py

Indentation matters in Python because it is used for defining scopes, which eliminates the use of braces {}. Please make sure the code at the same level are indented the same amount represented by tab or space characters.

Formatted Printing

You can print numbers using a format string similar to the printf function in C.

>>>a = 0.1
>>> a
0.1
>>>print "%20.19f" % a
0.1000000000000000056

This example show the limit of the floating point number representation. 1/10 can not be represented precisely because it is not a power of two. Even though the console prints a as 0.1 the underlying representation is almost 0.1.

Numerical Computation

>>> import math
>>> a = math.sqrt(2)
>>> a
1.4142135623730951
>>> a ** 2
2.0000000000000004

The following example will run forever till the result overflows the registers because x will never become exactly 1.0 because the representation of 0.1 is an approximation (with an error).

>>> x = 0.0
>>> while not x == 1.0:
...     x = x + 0.1 
...     print("x=%19.17g" % (x))
... 
x=0.10000000000000001
x=0.20000000000000001
x=0.30000000000000004
x=0.40000000000000002
x=                0.5
x=0.59999999999999998
x=0.69999999999999996

The moral of the lesson is: do not compare floating point numbers for strict equality. Alternatively we can calculate the distance between the two numbers and stop when the distance is short enough (less than a threshold).

>>> x = 0.0
>>> while abs(x - 1.0) > 1e-8:
...     x = x + 0.1
...     print("x=%19.17g" % (x))
... 
x=0.10000000000000001
x=0.20000000000000001
x=0.30000000000000004
x=0.40000000000000002
x=                0.5
x=0.59999999999999998
x=0.69999999999999996
x=0.79999999999999993
x=0.89999999999999991
x=0.99999999999999989

If we know exactly how many iterations to perform we can use a for loop as shown below. Note that range(1, 11) => [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].

>>> for i in range(1, 11):
...     x = i * 0.1
...     print("x=%19.17g" % (x))
... 
x=0.10000000000000001
x=0.20000000000000001
x=0.30000000000000004
x=0.40000000000000002
x=                0.5
x=0.60000000000000009
x=0.70000000000000007
x=0.80000000000000004
x=0.90000000000000002
x=                  1

What does the following piece of code compute? Will the loop ever finish? Why?

>>> eps = 1.0
>>> while 1.0 + eps > 1.0:
...     eps = eps / 2.0
... 
>>> print eps
1.11022302463e-16

The value eps has will always be a power of 2, therefore there will not be any rounding off errors. However, as the while loop continues eps gets smaller and smaller. Eventually eps become so small it becomes nothing comparing to 1.0, which means adding it to 1.0 will result in 1.0. Recall that the machine epsilon is the smallest value when added to 1.0 will result in value different than 1.0. If a number is less than the machine epsilon and it will cause no effect when it is added to 1.0. sys.float_info tells us that the machine epsilon is 2.220446049250313e-16, which is slightly less than 2x1.11022302463e-16=2.22044604926e-16. This is why when the loop ends eps equals 1.11022302463e-16.

Symbolic Computation

The SymPy library for Python allows us to manipulate variables symbolically (as symbols). The following code snippets demonstrates some basic functionality of SymPy.

>>> from sympy import Symbol
>>> x = Symbol('x');
>>> x+2*x+3*x
6*x
>>> y=x+2*x+3*x
>>> y
6*x
>>> y.subs(x, 2)
12
>>> from sympy import diff
>>> diff(y, x)
6
>>> from sympy import sin, cos
>>> diff(sin(x), x)
cos(x)
>>> diff(cos(x), x)
-sin(x)
>>> y = 4 - 1/x
>>> y.diff()
x**(-2)
>>> y.diff().subs(x, 0.5)
4.00000000000000
>>> from sympy import series
>>> sin(x).series(x, 0)
x - x**3/6 + x**5/120 + O(x**6)
>>> sin(x).series(x, 0, 10)
x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)

Numpy Array and Linear Algebra

Sources:

• Numpy Tutorial
• Scipy Tutorial

A array in numpy is multidimensional table of elements of the same type. The elements are indexed by a tuple of positive integers. The name of the class for numpy arrays is ndarray, a.k.a array. Each array object has a number of attributes:

• mdim: number of dimensions
• shape: a tuple of integers indicating the size of the array in each dimension. A matrix with m rows and n columns will have (m, n) as its shape.
• size: the total number elements in an array.
• dtype: the type of elements in an array.

You can install numpy and scipy on Ubuntu as follows:

sudo apt-get install python-numpy python-scipy python-matplotlib python-sympy

The elements in an array can be referenced using different syntax as shown in the examples below.

>>> from numpy import *
>>> a = arange(15).reshape(3, 5)
>>> a
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14]])
>>> type(a)
<type 'numpy.ndarray'>
>>> a.shape
(3, 5)
>>> (rows, cols) = a.shape
>>> rows
3
>>> cols
5
>>> a.ndim
2
>>> a.size
15
>>> a[2, 3]
13
>>> a[2][3]
13
>>> a[-1]
array([10, 11, 12, 13, 14])
>>> a[-2]
array([5, 6, 7, 8, 9])
>>> a[-2:]
array([[ 5,  6,  7,  8,  9],                                                 
       [10, 11, 12, 13, 14]])
>>> a[2:]
array([[10, 11, 12, 13, 14]])  
>>> a[:-3]
array([], shape=(0, 5), dtype=int64)
>>> a[:]
array([[ 0,  1,  2,  3,  4],                                                 
       [ 5,  6,  7,  8,  9],                                                 
       [10, 11, 12, 13, 14]]) 
>>> a[1, ...]
array([5, 6, 7, 8, 9]) 
>>> a[:, 0]
array([ 0,  5, 10])

A numpy array can be created from a regular Python list or tuple using the array function. The type of the array elements depends on that of the elements in the sequences.

>>> a = array([[1, 2], [2, 3]])
>>> a
array([[1, 2],
       [2, 3]])
>>> a = array(((4, 5), (6, 7)))
>>> a
array([[4, 5],
       [6, 7]])
>>> a.dtype
dtype('int64')
>>>>>> b = array([(1.2, 1.3), (1.4, 1.5)])
>>> b
array([[ 1.2,  1.3],
       [ 1.4,  1.5]])
>>> b.dtype
dtype('float64')

Numpy includes functions that create arrays with default values. The eye() function creates a identity matrix and the identity() performs a similar function. The copy function clones an array object includes the data it contains.

>>> zeros((2, 3))
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
>>> ones((3, 4))
array([[ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.]])
>>> empty((2, 3))
array([[  2.68156159e+154,   2.68156159e+154,   2.68156242e+154],
       [  2.68156159e+154,   2.68156159e+154,   2.68156159e+154]])
>>> eye(3, 3)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> identity(3, float)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

Arithmetic operations on numpy arrays are element-wise operations. The product * operator does a element-wise multiplication. The dot function performs matrix multiplications.

>>> a
array([[1, 2],
       [3, 4]])
>>> b = a.copy()
>>> b
array([[1, 2],
       [3, 4]])
>>> a*b
array([[ 1,  4],
       [ 9, 16]])
>>> dot(a, b)
array([[ 7, 10],
       [15, 22]])

Numpy also includes a number of basic linear algebra functions.

>>> from numpy.linalg import *
>>> a = array([[1.0, 2.0], [3.0, 4.0]])
>>> a
array([[ 1.,  2.],
       [ 3.,  4.]])
>>> a.transpose()
array([[ 1.,  3.],
       [ 2.,  4.]])
>>> inv(a)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])
>>> y = array([[5.], [7.]])
>>> solve(a, y)
array([[-3.],
       [ 4.]])

IPython Notebook

The IPython Notebook is an interactive environment that allows coding, executing, documentation of your Python commands.

Resources:

• How to share an IPython notebook
• A gallery of IPython notebooks
• A matplotlib tutorial
• Exploratory computing with Python