This IPython notebook is a direct translation of Chapter 1 of

- D. J. Higham and N.J. Higham.

MATLAB Guide, Second edition,

Society for Industrial and Applied Mathematics, Philadelphia, PA, USA,

2005, ISBN 0-89871-578-4

with MATLAB code converted to equivalent Python/Numpy (and IPython/Matplotlib/Scipy) code. There are a small number of additions in the text to address Python/Numpy/IPython differences and additions. You can find out more about IPython and the IPython Notebook here and in reference [2]. All the errors here are, of course, mine. I hope you find it useful.

Don MacMillen: don (sometimes) at macmillen dot net (see also this blog entry)

The best way to learn Python/Numpy is by trying it yourself, and hence we begin with a whirlwind tour. Working through the examples below will give you a feel for the way that Python/Numpy operates and an appreciation of its power and ﬂexibility.

The tutorial is entirely independent of the rest of the book—all the Python/Numpy features introduced are discussed in greater detail in the subsequent chapters. Indeed, in order to keep this chapter brief, we have not explained all the functions used here. You can use the index in [1] to ﬁnd out more about particular topics that interest you.

We will be using the IPython interactive shell as well as the IPython notebook (which is what you are reading now) throughout this tutorial. You can install these by following the directions at IPython install. This tutorial contains commands for you to type at the IPython command line. Alternatively, if you are viewing this notebook in a "live" mode served from a local IPython notebook server, you can directly modify any of the cells and rerun the cell. Just be aware that in some cases you may need to choose the "Run All" option from the cell drop down menu to satisfy intercell dependences.

In the last part of this brief tutorial we give examples of script files and functions. These ﬁles are short, so you can type them in quickly or cut and paste them into an IPython shell (using the %cpaste magic), or you can just modify values in place in the IPython notebook and re-evaluate the cell. You should experiment as you proceed, keeping the following points in mind.

- Upper and lower case characters are not equivalent (Python/Numpy are case sensitive).
- Typing the name of a variable will cause IPython to display its current value.
- Python/Numpy uses parentheses, ( ), square brackets, [ ], and curly braces, { }, and these are not interchangeable.
- In the IPython shell, the up arrow and down arrow keys can be used to scroll through your previous commands. Also, an old command can be recalled by typing the first few characters followed by up arrow. You can also use ctrl-p (type p while holding the control key down) to scroll back through the command history stack.
- You can type help(topic) to access online help on the command, function, or object topic. Note that hyperlinks, indicated by underlines, are provided that will take you to related help items and the Help browser.
- If you press the tab key after partially typing a function or variable name, IPython will attempt to complete it, offering you a selection of choices if there is more than one possible completion.
- You can quit IPython by typing exit or quit or ctrl-d twice.

Having entered the IPython command, you should work through this tutorial by typing in the text that appears after the IPython prompt, 'In [n]:' where n is a number which indicates the command's location in the command stack, in the command window. After showing you what to type, we display the output that is produced. We begin with arrays. In native Python arrays are lists denoted with the square bracket syntax. However, we want Numpy arrays, so first you import the numpy module and shorten the prefix to "np" (this is a common usage convention)

In [1]:

```
%matplotlib inline
import numpy as np
a = np.array([1, 2, 3])
a
```

Out[1]:

array([1, 2, 3])

If you type in the three lines listed in the box 'In[1]:' above into the IPython shell, then you will see the result listed in 'Out[1]:'

This example sets up a 3 element array, which is also called a vector. In some other languages, a distinction is made between a row vector and a column vector, but that is not the case in Numpy. If you define a new vector c

In [2]:

```
c = np.array([4, 5, 6])
c
```

Out[2]:

array([4, 5, 6])

Now you can multiply the arrays a and c:

In [3]:

```
a*c
```

Out[3]:

array([ 4, 10, 18])

In [4]:

```
_
```

Out[4]:

array([ 4, 10, 18])

In [5]:

```
np.dot(a, c)
```

Out[5]:

32

In [6]:

```
A = np.outer(c, a)
A
```

Out[6]:

array([[ 4, 8, 12], [ 5, 10, 15], [ 6, 12, 18]])

Here the answer is a 3-by-3 array that has been assigned to A.

The product a * a, since the '*' operation is element-wise, is equivalent to squaring a, which uses the '**' operator

In [7]:

```
a * a
```

Out[7]:

array([1, 4, 9])

In [8]:

```
b = a ** 2
b
```

Out[8]:

array([1, 4, 9])

Arithmetic operations on matrices and vectors come in two distinct forms. Array sense operations are deﬁned to act elementwise and, as mentioned before, are the default behavior of Numpy.

Matrix sense operations are based on the normal rules of linear algebra and are obtained with either matrix objects, or calling the specific matrix function from numpy. For matrix objects, the usual symbols +, -, *.

In [9]:

```
am = np.matrix(a)
cm = np.matrix(c)
am * cm.T
```

Out[9]:

matrix([[32]])

In [10]:

```
g = np.ones( (3, 3, 3) )
g
```

Out[10]:

array([[[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]], [[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]], [[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]]])

If you try to turn this into a matrix object, you will get an error

In [11]:

```
#h = np.matrix(g)
print("Error has a problem???")
```

Error has a problem???

Here you see how errors are displayed in Python. You see the call stack from the top to the bottom and the error message that is finally displayed.

Numpy has many mathematical functions that operate on arrays element wise when given an array argument. For example,

In [12]:

```
np.exp(a)
```

Out[12]:

array([ 2.71828183, 7.3890561 , 20.08553692])

In [13]:

```
np.log(_)
```

Out[13]:

array([ 1., 2., 3.])

In [14]:

```
np.sqrt(a)
```

Out[14]:

array([ 1. , 1.41421356, 1.73205081])

In [15]:

```
np.set_printoptions(precision=4)
print (np.sqrt(a))
```

[ 1. 1.4142 1.7321]

You can find out much more about the function set_printoptions by using IPython's help command

In [16]:

```
help(np.set_printoptions)
```

You set precision back to 8 to get the default behavior

In [17]:

```
np.set_printoptions(precision=8)
np.sqrt(a)
```

Out[17]:

array([ 1. , 1.41421356, 1.73205081])

In [18]:

```
2 ** -24
```

Out[18]:

5.960464477539063e-08

Various data analysis functions are also available.

In [19]:

```
b.sum()
```

Out[19]:

14

In [20]:

```
c.mean()
```

Out[20]:

5.0

In [21]:

```
dir(a)
```

Out[21]:

['T', '__abs__', '__add__', '__and__', '__array__', '__array_finalize__', '__array_interface__', '__array_prepare__', '__array_priority__', '__array_struct__', '__array_wrap__', '__bool__', '__class__', '__contains__', '__copy__', '__deepcopy__', '__delattr__', '__delitem__', '__dir__', '__divmod__', '__doc__', '__eq__', '__float__', '__floordiv__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__iadd__', '__iand__', '__ifloordiv__', '__ilshift__', '__imatmul__', '__imod__', '__imul__', '__index__', '__init__', '__int__', '__invert__', '__ior__', '__ipow__', '__irshift__', '__isub__', '__iter__', '__itruediv__', '__ixor__', '__le__', '__len__', '__lshift__', '__lt__', '__matmul__', '__mod__', '__mul__', '__ne__', '__neg__', '__new__', '__or__', '__pos__', '__pow__', '__radd__', '__rand__', '__rdivmod__', '__reduce__', '__reduce_ex__', '__repr__', '__rfloordiv__', '__rlshift__', '__rmatmul__', '__rmod__', '__rmul__', '__ror__', '__rpow__', '__rrshift__', '__rshift__', '__rsub__', '__rtruediv__', '__rxor__', '__setattr__', '__setitem__', '__setstate__', '__sizeof__', '__str__', '__sub__', '__subclasshook__', '__truediv__', '__xor__', 'all', 'any', 'argmax', 'argmin', 'argpartition', 'argsort', 'astype', 'base', 'byteswap', 'choose', 'clip', 'compress', 'conj', 'conjugate', 'copy', 'ctypes', 'cumprod', 'cumsum', 'data', 'diagonal', 'dot', 'dtype', 'dump', 'dumps', 'fill', 'flags', 'flat', 'flatten', 'getfield', 'imag', 'item', 'itemset', 'itemsize', 'max', 'mean', 'min', 'nbytes', 'ndim', 'newbyteorder', 'nonzero', 'partition', 'prod', 'ptp', 'put', 'ravel', 'real', 'repeat', 'reshape', 'resize', 'round', 'searchsorted', 'setfield', 'setflags', 'shape', 'size', 'sort', 'squeeze', 'std', 'strides', 'sum', 'swapaxes', 'take', 'tobytes', 'tofile', 'tolist', 'tostring', 'trace', 'transpose', 'var', 'view']

To find out more about any of these methods, use the help command

In [22]:

```
help(a.argmax)
```

In [23]:

```
np.pi
```

Out[23]:

3.141592653589793

In [24]:

```
_
```

Out[24]:

3.141592653589793

In [25]:

```
y = np.tan(np.pi/6)
y
```

Out[25]:

0.57735026918962562

The Variable np.pi is a permanent Variable with value $\pi$. The variable "_" always contains the most recent unassigned expression, as described earlier, so after the assignment to y, "_" still holds the Value $\pi$.

You may set up a two dimensional array by concatenating columns using the np.c_[ ] array notation

In [26]:

```
B = np.c_[[-3, 2, -1], [0, 5, 4], [-1, -7, 8]]
B
```

Out[26]:

array([[-3, 0, -1], [ 2, 5, -7], [-1, 4, 8]])

In [27]:

```
import numpy.linalg as nl
x = nl.solve(B, c)
x
```

Out[27]:

array([-1.29953917, 1.37788018, -0.10138249])

In [28]:

```
nl.norm(np.dot(B,x) - c) / (nl.norm(B) * nl.norm(x))
```

Out[28]:

7.2040771407550849e-17

Some times we see a 0 here, but we usually expect to see something nonzero because of rounding errors.

The eigenvalues of B can be found using eig from the numpy.linalg module

In [29]:

```
np.set_printoptions(precision=5) # make it the m*lab default
w, _ = nl.eig(B)
w
```

Out[29]:

array([-3.13605+0.j , 6.56803+5.10454j, 6.56803-5.10454j])

In [30]:

```
w, v = nl.eig(B)
v
```

Out[30]:

array([[ 0.98290+0.j , -0.03854-0.03928j, -0.03854+0.03928j], [-0.12656+0.j , -0.80053+0.j , -0.80053-0.j ], [ 0.13372+0.j , 0.16831+0.57254j, 0.16831-0.57254j]])

The columns of v are the eigenvectors of B.

To get vectors of evenly spaced values, use the np.arange function

In [31]:

```
v = np.arange(6)
v
```

Out[31]:

array([0, 1, 2, 3, 4, 5])

Note that the default starting value is zero. Also, as noted earlier, all Numpy arrays start at index zero and end at index n - 1 for an array of size n.

Nonunit increments can be specified by a third number, often called the 'stride'.

In [32]:

```
w = np.arange(2, 10, 3)
w
```

Out[32]:

array([2, 5, 8])

In [33]:

```
y = np.arange(1, 0, -0.25)
y
```

Out[33]:

array([ 1. , 0.75, 0.5 , 0.25])

*less* than 0 (but greater than -(0.25 + delta) ). Here we just use -0.001.

In [34]:

```
y = np.arange(1, -0.001, -0.25)
y
```

Out[34]:

array([ 1. , 0.75, 0.5 , 0.25, 0. ])

In [35]:

```
C = np.c_[A, [8, 9, 10]]
C
```

Out[35]:

array([[ 4, 8, 12, 8], [ 5, 10, 15, 9], [ 6, 12, 18, 10]])

In [36]:

```
D = np.r_[B, B, B]
D
```

Out[36]:

array([[-3, 0, -1], [ 2, 5, -7], [-1, 4, 8], [-3, 0, -1], [ 2, 5, -7], [-1, 4, 8], [-3, 0, -1], [ 2, 5, -7], [-1, 4, 8]])

In [37]:

```
E = np.vstack((B, B, a))
E
```

Out[37]:

array([[-3, 0, -1], [ 2, 5, -7], [-1, 4, 8], [-3, 0, -1], [ 2, 5, -7], [-1, 4, 8], [ 1, 2, 3]])

Notice well that np.vstack takes a *single* argument that is a tuple of arrays. That way you can stack as many arrays as needed with a single call to vstack or just use the r_[ ] notation.

The element in row i and column j of the matrix C (where i and j always start at 0) can be accessed as c[i, j]

In [38]:

```
C[1, 2]
```

Out[38]:

15

*slicing*.

In [39]:

```
C[1:3, 0:2]
```

Out[39]:

array([[ 5, 10], [ 6, 12]])

In [40]:

```
I3 = np.eye(3)
I3
```

Out[40]:

array([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]])

In [41]:

```
Y = np.zeros((3,5))
Y
```

Out[41]:

array([[ 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0.]])

In [42]:

```
Z = np.ones((2,2))
Z
```

Out[42]:

array([[ 1., 1.], [ 1., 1.]])

Note that these functions take a single argument which is a tuple (except for np.eye, which has a different interface). The first argument of the tuple speciﬁes the size of the first dimension of the array and so on.

The methods rand and randm on an object returned from numpy.random.RandomState() do not take tuples. They are convenience functions that have a similar interface to another math envirionment. These two methods generate random entries from the uniform distribution over [0,1] and the normal (0,1) (ie zero mean with unit variance) distribution, respectively.

Note that it is always good for you to intialize a new instance of RandomState when using random numbers. This ensures that the random stream is for you alone, and you can reset the seed to make the sequence repeatable. Here we set the starting seed to 20.

In [43]:

```
rn = np.random.RandomState() # initialize a new RandomState object
rn.seed(20)
F = rn.rand(3, 3)
F
```

Out[43]:

array([[ 0.58813, 0.89771, 0.89153], [ 0.81584, 0.03589, 0.69176], [ 0.37868, 0.51851, 0.65795]])

In [44]:

```
G = rn.randn(1, 5)
G
```

Out[44]:

array([[-0.62064, -0.83453, 0.91636, 0.70784, 0.41968]])

In [45]:

```
%who
```

A B C D E F G I3 Y Z a am b c cm g nl rn v w x y

The %whos magic will additionally give the types and some additional information

In [46]:

```
%whos
```

which approximates the golden ratio, (1 + $\sqrt{5}$)/2. The evaluation is done from the bottom up:

In [47]:

```
g = 2.
for k in range(10):
g = 1 + 1 / g
g
```

Out[47]:

1.6180555555555556

In [48]:

```
import scipy.constants as sconst
sconst.golden
```

Out[48]:

1.618033988749895

Loops involving the while statement can be found later in this tutorial.

The plot function from the matplotlib module pylab produce two dimensional pictures. You can assess these ploting and graphics functions by importing this module. It is a common practice to rename it to plt, as in the following example.

In [49]:

```
#import mpld3
#mpld3.enable_notebook()
import matplotlib.pylab as plt
t = np.arange(0, 1.005, 0.005)
z = np.exp(10 * t * (t - 1)) * np.sin(12 * np.pi * t)
plt.plot(t, z)
plt.title("Figure 1.2. Basic 2D picture produced by plt.plot")
```

Out[49]:

<matplotlib.text.Text at 0x111f56c88>

Here, plt.plot(t, z) joins the points t[i], z[i] using the default solid linetype. Matplotlib opens a figure window in which the picture is displayed. In this IPython notebook, the matplotlib figures are included 'in line' since the notebook was invoked with the "ipython notebook --pylab=inline" (the --pylab=inline option will be deprecated in the 3.0 notebook in favor of the **%matplotlib inline** cell magic) command. When you are using IPython as a command line interpreter, after opening a plot figure window, you can close it in the normal way, ie by clicking on the x in the window title bar.

You can produce a histogram with the function plt.hist

In [50]:

```
plt.hist(np.random.randn(1000))
plt.title("Figure 1.3. Histogram produced by plt.hist.")
plt.show()
```

Here, hist is given 1000 points from the normal (0, 1) random number generator

You are now ready for more challenging computations. A random Fibonacci sequency {$x_{n}$} is generated by choosing $x_{1}$ and $x_{2}$ and setting

$$ x_{n + 1} = x_{n} \pm x_{n -1}, \; n \ge 2 $$Here, the $\pm$ indicates that + and - must have equal probability of being chosen. Viswanath [121] listed in [1] analyzed this recurrence and showed that, with probability 1, for large n the quantity $|x_{n}|$ increases like a multiple of $c^{n}$, where c = 1.13198824... (see also [25]). You can test Viswanath's result as follows:

In [51]:

```
%reset -f
import numpy as np
import pylab as plt
rd = np.random.RandomState()
rd.seed(100)
x = np.zeros(1000)
x[0:2] = 1, 2
for n in range(1, 999):
x[n + 1] = x[n] + np.sign(rd.rand(1) - 0.5) * x[n - 1]
xx = np.arange(1, 1001)
plt.semilogy(xx, np.abs(x))
c = 1.13198824
plt.hold(True)
plt.semilogy(xx, c ** xx)
plt.title("Figure 1.4. Growth of a random Fibonacci sequence")
plt.hold(False)
```

Here, %reset -f removed all variables and imported modules from the workspace, which is why we needed to import numpy and pylab again. The for loop stores a random Fibonacci sequence in the array x; we preallocate x to the size needed and initial to zero using the np.zeros() function. The plt.semilogy function then plots n on the x-axis against $|X|$ on the y-axis, with logarithmic scaling for the y-axis. Typing np.hold(True) tells matplotlib to superimpose the next picture on top of the current one. The second semilogy plot produces a line of slope c. The overall picture, shown in Figure 1.4, is consistent with Viswanath’s theory.

We can make the above code into a command by writing it out to a file. If you cut and paste it into a file called fib.py, you can then run it in IPython by typing "run fib" or "run fib.py" to reproduce the above graph.

However, you can experiment directly in this IPython notebook by changing any of the values and then hitting the 'play' button above.

Our next example involved the Collatz iteration, which, given a positive integer $x_{1}$, has the form $x_{k + 1} = f(x_{k})$, where

$$f(x) = \left\{ \begin{array}{lr} 3x + 1 & : x \: \% \: 2 == 0\\ x/2 & : x \: \% \: 2 == 1 \end{array} \right. $$Here $x \: \% \: y$ is the modulus function of x and y, sometimes also written as mod(x, y). It is the remainder of x when divided by y. In this case $x \: \% \: 2 == 0$, if true, means that x is even and if $x \: \% \: 2 == 1$ is true, then x is odd. Note that these are perfectly good Python functions

In [52]:

```
x = 3
if x % 2 == 0:
print ('x is even')
else:
print ('x is odd')
```

x is odd

In [53]:

```
#COLLATZ Collatz iteration.
#n = int(raw_input(’Enter an integer bigger than 2: ’))
init_n = n = 27
narray = np.zeros(1000) # only a maximum of 1000 iterations
narray[0] = n
count = 0
while n != 1:
if n % 2 == 1: # Remainder modulo 2.
n = 3 * n + 1
else:
n = n / 2
count += 1
narray[count] = n # Store the current iterate
# Plot with * marker and solid line style.
# Only plot the non zero entries in narray
plt.plot(narray[narray != 0], '*-')
plt.title('Figure 1.5. Collatz iteration starting at %d' % init_n)
```

Out[53]:

<matplotlib.text.Text at 0x11a0887f0>

A couple of things you should note about the code above. This first is that it will fail if the number of iterations goes beyond 1000. For most of the inputs I have tried, it 'converged' in under 200 iterations. But you should try to break it! The second thing to note is that we only want to plot the non zero values. All the zero values at the end of the array add no new information and we want to drop them. That is easily accomplished by using the narray[narray != 0] syntax which selects only the non-zero values out of narray.

To investigate the Collatz problem further, the script collbar in the next listing plots a bar graph of the number of iterations required to reach the value 1, for starting Values 1,2,. . . ,29. The result is shown in Figure 1.6. For this picture, the function plt.grid adds grid lines that extend from the axis tick marks, while plt.title, plt.xlabel, and plt.ylabel add further information.

In [54]:

```
#COLLBAR Collatz iteration bar graph.
N = 29 # Use starting values 1,2,...,N.
niter = np.zeros(N); # Preallocate array.
for i in range(N):
count = 0
n = i + 1
while n != 1:
if n % 2 == 1:
n = 3 * n + 1
else:
n = n / 2
count += 1
niter[i] = count
left = np.arange(29)
plt.bar(left, niter) # Bar graph.
plt.grid() # Add horizontal and vertical grid lines.
plt.title('Figure 1.6. Col1atz iteration counts')
plt.xlabel('Starting value') # Label X-axis.
plt.ylabel('Number of iterations') #Label y-axis.
```

Out[54]:

<matplotlib.text.Text at 0x11a0889b0>

In [55]:

```
#MANDEL Mandelbrot set.
x = np.linspace(-2.1, 0.6, 301)
y = np.linspace(-1.1, 1.1, 301)
[X,Y] = np.meshgrid(x, y)
C = X + 1j * Y
Z_max = 1e6
it_max = 50
Z = C
for k in range(it_max):
Z = Z ** 2 + C
plt.contourf(x, y, np.float64(np.abs(Z) < Z_max))
plt.title("Figure 1.7 Mandelbrot Set")
plt.show()
```

Next we solve the ordinary differential equation (ODE) system

$$ d/dt \; y_{1}(t) = 10 (y_{2}(t) - y_{1}(t)),$$$$ $$$$ d/dt \; y_{2}(t) = 28 y_{1}(t) - y_{2}(t) - y_{1}(t) y_{3}(t),$$$$ $$$$ d/dt \; y_{3}(t) = y_{1}(t) y_{2}(t) - 8 y_{3}(t) / 3. $$This is an example from the Lorenz equations family; see [H1] listed in [1]. We take initial conditions y(0) = [0, 1, 0] and solve over $0 \le t \le 50$. The next listing, titled Lorentz, is an example of a Python function. Given t and y, this function returns the right-hand side of the ODE as the vector **yprime**. This is the form required by Scipy’s ODE solving functions. The rest of the listing uses the scipy function odeint to solve the ODE numerically and then produces the (y1, y3) phase plane plot shown in Figure 1.8. You can try different values of the constants defining the derivatives to see what happens. For instance, change the 3 to an 8 in lorenzde and observe how the plot changes

In [56]:

```
import scipy.integrate as si
def lorenzde(y, t):
'''LORENZDE Lorenz equations.
YPRIME = LORENZDE(Y,T)
'''
yprime = np.array([10. * (y[1] - y[0]), 28. * y[0] - y[1] - y[0] * y[2],
y[0] * y[1] - 8. * y[2] / 3.])
return yprime
#lrun ODE solving example: Lortez.
t = np.arange(0, 50.01, .01) # time points on which to solve
yzero = np.array([0., 1., 0.])
print (len(yzero))
y = si.odeint(lorenzde, yzero, t)
plt.plot(y[:, 0], y[:, 2])
plt.xlabel('y_0')
plt.ylabel('y_2')
plt.title('Figure 1.8 Lorenz equations')
```

3

Out[56]:

<matplotlib.text.Text at 0x11c216240>