Introduction¶

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)

Chapter 1 A Brief Tutorial¶

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 flexibility.

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 find 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 files 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)

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

Now you can multiply the arrays a and c:

Here, you performed an element by element mulitpy of the two vectors, not the scalar or dot product that some other languages would have done. Notice that in this example, the product a*c was not explicitly assigned to any variable. When there is no explicit assignment, IPython automatically assigns the result of the expression to the variable named '_' (underscore) as you can see by the following:

If you wanted to get the dot product of the two vectors, use the dot function from numpy, as in the following

Inputs to Python functions are specified after the function name and within parentheses, as you just saw above with the 'dot' function from the numpy module. You may also form the outer product of the two vectors by calling another function from numpy

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

Arithmetic operations on matrices and vectors come in two distinct forms. Array sense operations are defined 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 +, -, *.

Notice that you had to take the transpose of c in order to make it into a column vector. For numpy matrix objects, there is a difference between row and column vectors. Also, matrix object can only have one or two dimensions. For instance, the numpy function 'ones' takes a tuple and returns a numpy array of that dimensionality. The following will give us a 3 by 3 by 3 array of ones.

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

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,

By default, Numpy displays floating point numbers to 8 decimal digits, but always stores numbers and computes to the equivalent of 16 decimal digits. The output format can be changed using the set_printoptions function

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

You set precision back to 8 to get the default behavior

Large or small numbers are displayed in exponential notation, with a power of 10 scale factor preceded by e:

Various data analysis functions are also available.

Here you see an example of the object oriented nature of Numpy. All Numpy arrays are of type "ndarray" and have many methods associated with them. You can list them all by using a 'dir' command in IPython, use the following

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

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

At the heart of Numpy is a powerful range of linear algebra functions. For example, recalling that c is a 3-by-1 Vector, you may wish to solve the linear system $B * x = c$. This can be done with the solve function from the numpy.linalg module

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

Notice two things about this code snippet. First, nl.eig returns two values, w are the eigenvalues and v would be the eigenvectors. Since in this instance we did not want the eigenvectors we set it to the underscore '_', which is a common idiom in Python. Finally, the solution shows the j is the imaginary unit, $\sqrt{-1}$. To get the eigenvectors you can do the calculation again

The columns of v are the eigenvectors of B.

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

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'.

This last example illustrates a peculiarity of Python and numpy, that the ending value or index is actually one position short of what you might expect. In the previous example, to have the array actually increment down to zero, you have to specify an end value that is strictly less than 0 (but greater than -(0.25 + delta) ). Here we just use -0.001.

You may construct big matrices out of little ones by using the numpy functions hstack and vstack and the previously mentioned c_[ ] and r_[ ] notations

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]

More generically, C[i1:i2, j1:j2] picks out the submatrix formed by the intersection of rows i1 to i2-1 and columns j1 to j2-1. This type of operation on arrays is generally called slicing.

You can build certain types of matrices automatically. For example, identities and matrices of zeros and ones can be constructed with eye, zeros, and ones:

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 specifies 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.

At this point several variables have been created in the workspace. You can obtain a list with the %who command, where the % indicates an IPython 'magic' command.

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

Like most languages Python has loop contructs. The following example uses a for loop to evaluate the continued fraction $\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{1 + \cfrac{1}{1 + \cfrac{1}{1 + \cfrac{1}{1}}}}}}}}}}$

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

Many constants can be found in the scipy submodule constants. You can access them in the following way

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.

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

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:

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

But now returning to our quest, the equation above in words means: if $x$ is odd, replace it by $3x + 1$, and if $x$ is even, halve it. It has been conjectured that this iteration will always lead to a value of 1 (and hence thereafter cycle between 4, 2, and 1) whatever starting value $x_{1}$ is chosen. There is ample computational evidence to support this conjecture, which is variously known as the Collatz problem, the $3x + 1$ problem, the Syracuse problem, Kakutani’s problem, Hasse’s algorithm, and Ulam’s problem. However, a rigorous proof has so far eluded mathematicians. For further details, see [63] listed in [1] or type “Collatz problem” into your favorite Web search engine. You can investigate the conjecture by creating the python script file collatz.py shown below. In this file a while loop and an if statement are used to implement the iteration. You can run this in the notebook simply by changing the value of n and hitting the play button. Or, you can save this to a file named collatz.py (remembering to uncommment out the line '#n = int(raw_input..." and commenting out the line "n = 27". Then, from the IPython shell prompt, type "run collatz". The input command prompts you for a starting value. The appropriate response is to type an integer and then hit return.

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.

The Well-known and much studied Mandelbrot set can be approximated graphically in just a few lines of Python/Numpy. It is defined as the set of points c in the complex plane for which the sequence generated by the map $z \mapsto z^{2} + c$, starting with $z = c$, remains bounded [91, Chap. 14] listed in [1]. The script mandel in the next listing produces the plot of the Mandelbrot set shown in Figure 1.7. The script contains calls to np.linspace of the form np.linspace (a, b ,n), which generate an equally spaced vector of n values between a and b. The meshgrid and complex functions are used to construct a matrix C that represents the rectangular region of interest in the complex plane. The plot itself is produced by plt.contourf, which plots a filled contour. The expression $abs(Z)<Z\_max$ in the call to contourf detects points that have not exceeded the threshold Z_max and that are therefore assumed to lie in the Mandelbrot set; the double function is applied in order to convert the resulting logical array to numeric form. You can experiment with mandel by changing the region that is plotted, via the linspace calls, the number of iterations it_max, and the threshold Z_max. Also note that on some runs we may see a RuntimeWarning printed becase we have an overflow on z. This can be silenced with a context manager but we will wait to introduce detailed control of errors and warnings.

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

Now we give an example of a recursive function, that is, a function that calls itself. The Sierpinski gasket [90, Sec. 2.2] listed in [1] is based on the following process. Given a triangle with vertices $P_{a}$, $P_{b}$, and $P_{c}$, We remove the triangle with vertices at the midpoints of the edges, $(P_{a} + P_{b})/2$, $(P_{b} + P_{c}) / 2$, and $(P_{c} + P_{a}) / 2$. This removes the “middle quarter” of the triangle, as illustrated in Figure 1.9. (The code in the function 'gasket' will be explained below)

Effectively, we have replaced the original triangle with three “subtriangles”. We can now apply the middle quarter removal process to each of these subtriangles to generate nine subsubtriangles, and so on. The Sierpinski gasket is the set of all points that are never removed by repeated application of this process. The function gasket in the listing above implements the removal process. The input arguments $P\_{a}$, $P\_{b}$, and $P\_{c}$ define the vertices of the triangle and level specifies how many times the process is to be applied. If level is nonzero then gasket calls itself three times with level reduced by 1, once for each of the three subtriangles. When level finally reaches zero, the recursion 'bottoms-out' and the appropriate triangle is drawn. The following code generates Figure 1.10.

(Figure 1.9 was generated in the same way with level = 1.) In the last line, the call to axis makes the units of the $x-$ and $y—$axes equal and turns off the axes and their labels. You should experiment with different initial Vertices $P_{a}$, $P_{b}$, and $P_{c}$, and diiferent levels of recursion, but keep in mind that setting level bigger than 8 may overstretch either your patience or your computer’s resources.

The Sierpinski gasket can also be generated by playing Barnsley’s “chaos game” [90, Sec. 1.3] listed in [1]. We choose one of the vertices of a triangle as a starting point. Then we pick one of the three vertices at random, take the midpoint of the line joining this vertex with the starting point and plot this new point. Then we take the midpoint of the line joining this point and a randomly chosen vertex as the next point, which is plotted, and the process continues. The script barnsley in Listing 1.8 implements the game. Figure 1.11 shows the result of choosing 1000 iterations:

Try experimenting with the number of points, $n$, the type and size of marker in the plot command, and the location of the starting point.

We finish with the listing of sweep, which generates a volume—swept three—dimensional (3D) object; see Figure 1.12. Here, the command surf (X,Y,Z) creates a 3D surface where the height Z[i, j] is specified at the point (X[i, j] ,Y[i, j]) in the x-y plane. The script is not written in the most obvious fashion, which would use two nested for loops. Instead it is vectorized. To understand how it works you will need to be familiar with Chapter 5 and 21.4 of reference [1]. You can experiment with the script by changing the parameter N and the function that determines the variable radius: try replacing sqrt by other functions, such as log, sin, or abs.

[1] 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

[2] Fernando Pérez, Brian E. Granger, IPython: A System for Interactive Scientific Computing, Computing in Science and Engineering, vol. 9, no. 3, pp. 21-29, May/June 2007, doi:10.1109/MCSE.2007.53. URL: http://ipython.org