# 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:

``` Python
>>>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:

``` python
>>>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.

``` Python
>>> 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().

``` Python
>>> 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.

``` Python
>>> 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.

``` Python
>>> 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:

``` bash
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.

``` Python
>>>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

``` Python
>>> 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).

``` Python
>>> 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).

``` Python
>>> 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\].

``` Python
>>> 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?

``` Python
>>> 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.

``` python
>>> 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:

``` shell
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.

``` python
>>> 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.

``` python
>>> 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.

``` python
>>> 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.

``` python
>>> 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.

``` python
>>> 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

