A lot of computational algorithms are expressed using Linear Algebra terminology - vectors and matrices. This is thanks to the wide range of methods within Linear Algebra for solving the sort of problems that computers are good at solving!

Within Python, our first thought may be to represent a vector as a list. But there is a downside: lists do not naturally behave as vectors. For example:

Similarly, we cannot apply algebraic operations or functions to lists in a straightforward manner that matches our expectations.

However, there is a Python package - numpy - that does give us the behaviour we want.

numpy

numpy is used so frequently that in a lot of cases and online explanations you will see it abbreviated, using import numpy as np. Here we try to avoid that - auto-completion inside spyder means that the additional typing is trivial, and using the full name is clearer.

numpy defines a special type, an array, which can represent vectors, matrices, and other higher-rank objects. Unlike standard Python lists, an array can only contain objects of a single type. The notation to create these objects is straightforward: one easy way is to start with a list:

We see that the array objects behave as we would expect, and accessing elements is exactly the same as for a list. We can also perform other mathematical operations on the whole vector:

Think about these carefully.

1. The first case is straightforward: all elements of the vector are multiplied by a constant.
2. The second case applies a function to each element separately. numpy implements a version of most interesting mathematical functions, which are applied directly to each element.
3. The third case is elementwise multiplication of the vectors. The first component of the answer is the product of the first component of x_numpy with the first component of y_numpy. The second component of the answer is the product of the second component of x_numpy with the second component of y_numpy. We cannot use the * operator to represent matrix multiplication, but must use a function (see below; note that there will be an operator in Python 3.5+, but using it will mean your code is, for now, not widely useable).
4. The fourth case shows a combination of cases above. The answer is given by elementwise subtraction of 1, then squaring (elementwise) that result.

Defining a matrix can be done by applying the array function to a list of lists:

We see that for higher rank objects such as matrices, the operations are still performed elementwise.

To multiply a matrix by a vector, or a vector by a vector, in the standard linear algebra sense, we use the numpy.dot function:

Note that we have appeared to multiply a matrix by a row vector, and get a row vector back. This is because numpy does not distinguish between row and column vectors, so everything appears as a row vector. (You could define a $n \times 1$ array instead, but there is no advantage).

To actually check the shape and size of numpy arrays, you can directly check their attributes:

numpy contains a number of very efficient functions for working with arrays, for finding extreme values, and performing linear algebra tasks. Particular functions that are worth knowing, or starting from, are

• arange: constructs an array containing increasing integers
• linspace: constructs a linearly spaced array
• zeros and ones: constructs arrays containing just ones or zeros
• diag: extracts the diagonal of a matrix, or build a matrix with just diagonal entries
• mgrid: constructs matrices from vectors for 3d plots
• random.rand: constructs an array of random numbers.

Linear algebra

numpy also defines a number of linear algebra functions. However, a more comprehensive set of functions, which is better maintained and often more efficient, is given by scipy:

In addition to solving linear systems and computing determinants, you can also factorize matrices and generally do most linear algebra operations that you need. The scipy documentation is comprehensive, and has a specific section on Linear Algebra, as well as a section in the tutorial. Johansson also has a tutorial on scipy in general.

Working with files

Often we will want to work with data - constants, parameters, initial conditions, measurements, and so on. numpy provides ways to work with data stored in files - either reading them in or writing them out. A list of "File I/O routines" is available, but the two key routines are loadtxt and savetxt.

As a simple example we take our matrix A_numpy above and save it to a file:

We can then check the contents of that file (you should open the file on your machine to check):

Finally, we can read the contents of that file into a new variable and check that it matches:

Plotting

There are many Python plotting libraries depending on your purpose. However, the standard general-purpose library is matplotlib. This is often used through its pyplot interface.

This is a quick recap of the basic plotting commands, but using numpy as well.

This plotting interface is straightforward, but the results are not particularly nice. The following commands illustrate some of the ways of improving the plot:

Whilst most of the commands are self-explanatory, a brief note should be made of the strings line r'$x$'. These strings are in LaTeX format, which is the standard typesetting method for professional-level mathematics. The $ symbols surround mathematics. The r before the definition of the string says that the following string will be "raw": that backslash characters should be left alone. Then, special LaTeX commands have a backslash in front of them: here we use \pi and \sin. We can also use ^ to denote superscripts (used here), _ to denote subscripts, and use {} to group terms.

By combining these basic commands with other plotting types (semilogx and loglog, for example), most simple plots can be produced quickly.

Saving figures

If you want to save the figure to a file, instead of printing it to the screen, use the savefig command instead of the show command. For example, try:

We can then check the file on disk (you should open the file on your machine to check):

The type of the file is taken from the extension. Here we have used a png file, but svg and pdf output will also work.

Object-based approach

To get a more detailed control over the plot it's better to look at the objects that matplotlib is producing. Remember, when we talked about classes we said that it is an object with attributes and methods (functions) that are accessed using dot notation. Here are steps to completely control the plot.

First we define a figure object. We do not have to define the figure class - it is defined within matplotlib itself, along with a lot of useful methods. We call the constructor of the figure object in the same way as in the previous section, by calling pyplot.figure(). We can and will control its size (the units default to inches) by passing additional arguments to the constructor:

We will then define two axes on this figure. The numbers refers to the positions of the edges of the axes with respect to the figure window (between 0 and 1):

We will then add data to the both axes:

We will then set the range of the second axis:

Finally, we'll see what it looks like:

Each axis contains additional objects that can be modified.

Adding multiple axes by hand is often annoying (although sometimes necessary). There are a number of tools that can be used to simplify this in standard cases: add_subplot is the standard one. When you want a figure containing multiple subplots all the same size, with r rows and c columns, the command is add_subplot(r, c, <subplot_number>). For example:

The tight_layout function call at the end ensures that the axis labels and titles do not overlap with other subplots.

Higher dimensions

To plot three-dimensional objects, we need to modify the axis so that it knows a third dimension is required. To do this, we import another module and modify the command that sets up the axis object.

We can then construct, for example, a parametric spiral:

If we want to plot a surface, then we need to construct 2d arrays containing the locations of the $x$ and $y$ coordinates, and a 2d array containing the "height" of the surface. For structured data (ie, where the $x$ and $y$ coordinates lie on a regular grid) the meshgrid function helps. For example, the function

\begin{equation} \phi(x, y) = \sin^2 ( \pi x y ) \cos( 2 \pi y^2 ), \qquad x \in [0, 1], \quad y \in [0, 1], \end{equation}

would be plotted using

There are a lot of options to modify the appearance of this plot. Important ones include the colormap (note the US spelling), which requires importing the cm module from matplotlib, and the stride parameters changing the appearance of the grid. For example

Further reading

As noted earlier, the matplotlib documentation contains a lot of details, and the gallery contains a lot of examples that can be adapted to fit. There is also an extremely useful document as part of Johansson's lectures on scientific Python.

scipy

scipy is a package for scientific Python, and contains many functions that are essential for mathematics. It works particularly well with numpy. We briefly introduced it above for tackling Linear Algebra problems, but it also includes

• Scientific constants
• Integration and ODE solvers
• Interpolation
• Optimization and root finding
• Statistical functions

and much more.

Integration

The numerical quadrature problem involves solving the definite integral

\begin{equation} \int_a^b f(x) \, \text{d} x, \end{equation}

or a suitable generalization. scipy has a module, scipy.integrate, that includes a number of functions to solve these types of problems. For example, to solve

\begin{equation} I = \int_0^{\pi} \sin^2(x) \, \text{d} x, \end{equation}

the quad function can be used as:

The steps we have taken are:

1. Define the integrand by defining a function. This function takes the points at which the integrand is evaluated. By using numpy we can do this with a single command.
2. Import the quad function.
3. Call the quad function, passing the function defining the integrand, and the lower and upper limits.

The result we get back, as seen from the screen output, is not just $I$. It is a tuple containing both $I$, and also the accuracy with which quad believes it has computed the result. The quadrature is a numerical approximation, so can never be perfect. You should check this error estimate to ensure the result is "good enough" for your purposes.

We can also pass additional parameters if needed. Consider the problem

\begin{equation} I_a = \int_0^{\pi} \sin^2 (a x) \, \text{d} x. \end{equation}

If we wanted to solve this for many values of $a$, say $a = 1, 2, \dots, 5$, we could create a function taking a parameter, and then pass that parameter through:

Note that when passing the parameters using the args keyword argument, we put the parameters in a tuple. This shows how to pass more than one parameter: keep adding parameters to the argument list, and add them to the tuple. For example, to solve

\begin{equation} I_{a,b} = \int_0^{\pi} \sin^2(ax + b) \, \text{d} x \end{equation}

we write

Solving ODEs

There is a link between the solution of integrals and the solution of differential equations. Unfortunately, the numerical solution of an ODE is more complex than the solution of an integral. Fortunately, scipy contains a number of methods for these as well.

The methods in scipy solve ODEs of the form \begin{equation} \frac{\text{d} \vec{y}}{\text{d} t} = \vec{f} \left( \vec{y}, t \right), \qquad \vec{y}(0) = \vec{y}_0. \end{equation} For example, the ODE \begin{equation} \frac{\text{d} y}{\text{d} t} = e^{-t} - y, \qquad y(0) = 1 \end{equation} has $f(y, t) = e^{-t} - y$.

The method for using scipy is similar to the integration case.

1. Define a function that specifies the system, by defining the RHS.
2. Import the function that solves ODEs (odeint)
3. Call the function, passing the RHS function, the initial data $\vec{y}_0$, the times at which the solution is needed, and any parameters.

To solve our example, we use:

Note that the result for $y$ is not a vector, but a two dimensional array. This is because scipy will solve a general system of ODEs. This scalar case is a system of size $1$, but it still returns an array. To solve a system, the RHS function must take a vector for y, return a vector for dydt, and the initial data y0 must be a vector. All these vectors must be the same size.

The output is the numerical approximation to $y$ at the input times $t$, and can be immediately plotted:

Passing parameters is also similar to the integration case. For example, consider the problem

\begin{equation} \frac{\text{d}}{\text{d} t} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} -y + \alpha \\ x \end{pmatrix}, \qquad \begin{pmatrix} x \\ y \end{pmatrix}(0) = \begin{pmatrix} 1 \\ 0 \end{pmatrix}. \end{equation}

If $\alpha$ is zero, the solution is a circle in the $x, y$ plane. We solve this using odeint, denoting the state vector $\vec{z} = (x, y)^T$:

Further reading

Earlier we introduced scipy for Linear Algebra, and gave links there. Most of those links cover the full scipy package. The scipy documentation is comprehensive. Johansson also has a tutorial on scipy.

Exercise: Lorenz attractor

The Lorenz system is a set of ordinary differential equations which can be written

\begin{equation} \frac{\text{d} \vec{v}}{\text{d} \vec{t}} = \vec{f}(\vec{v}) \end{equation}

where the variables in the state vector $\vec{v}$ are

\begin{equation} \vec{v} = \begin{pmatrix} x(t) \\ y(t) \\ z(t) \end{pmatrix} \end{equation}

and the function defining the ODE is

\begin{equation} \vec{f} = \begin{pmatrix} \sigma \left( y(t) - x(t) \right) \\ x(t) \left( \rho - z(t) \right) - y(t) \\ x(t) y(t) - \beta z(t) \end{pmatrix}. \end{equation}

The parameters $\sigma, \rho, \beta$ are all real numbers.

Exercise 1

Write a function dvdt(v, t, params) that returns $\vec{f}$ given $\vec{v}, t$ and the parameters $\sigma, \rho, \beta$.

Exercise 2

Fix $\sigma=10, \beta=8/3$. Set initial data to be $\vec{v}(0) = \vec{1}$. Using scipy, specifically the odeint function of scipy.integrate, solve the Lorenz system up to $t=100$ for $\rho=13, 14, 15$ and $28$.

Plot your results in 3d, plotting $x, y, z$.

Exercise 3

Fix $\rho = 28$. Solve the Lorenz system twice, up to $t=40$, using the two different initial conditions $\vec{v}(0) = \vec{1}$ and $\vec{v}(0) = \vec{1} + \vec{10^{-5}}$.

Show four plots. Each plot should show the two solutions on the same axes, plotting $x, y$ and $z$. Each plot should show $10$ units of time, ie the first shows $t \in [0, 10]$, the second shows $t \in [10, 20]$, and so on.

This shows the sensitive dependence on initial conditions that is characteristic of chaotic behaviour.

Exercise: The shortest published Mathematical paper

A candidate for the shortest mathematical paper ever shows the following result:

\begin{equation} 27^5 + 84^5 + 110^5 + 133^5 = 144^5. \end{equation}

This is interesting as

This is a counterexample to a conjecture by Euler ... that at least $n$ $n$th powers are required to sum to an $n$th power, $n > 2$.

Exercise 1

Using Python, check the equation above is true.

Exercise 2

The more interesting statement in the paper is that

\begin{equation} 27^5 + 84^5 + 110^5 + 133^5 = 144^5. \end{equation}

[is] the smallest instance in which four fifth powers sum to a fifth power.

Interpreting "the smallest instance" to mean the solution where the right hand side term (the largest integer) is the smallest, we want to use Python to check this statement.

We are going to need to generate all possible combinations of four integers $a, b, c, d$ and test if $a^5 + b^5 + c^5 + d^5$ matches $e^5$ where $e$ is another integer.

The problem is the number of combinations grows very fast - the standard formula says that for a list of length $n$ there are

\begin{equation} \begin{pmatrix} n \\ k \end{pmatrix} = \frac{n!}{k! (n-k)!} \end{equation}

combinations of length $k$. For $k=4$ as needed here we will have $n (n-1) (n-2) (n-3) / 24$ combinations.

Show, by getting Python to compute the number of combinations $N = \begin{pmatrix} n \\ 4 \end{pmatrix}$ that $N$ grows roughly as $n^4$. To do this, plot the number of combinations and $n^4$ on a log-log scale. Restrict to $n \le 50$.

You may find the combinations function from the itertools package useful.

Exercise 3

With 17 million combinations to work with, we'll need to be a little careful how we compute. To check the interesting statement in the paper,

1. Construct a numpy array containing all integers in $1, \dots, 144$ to the fifth power.
2. Construct a list of all combinations of four elements from this array.
3. Construct a list of sums of all these combinations.
4. Loop over one list and check if the entry appears in the other list (ie, use the in keyword).

By printing out any entries that pass this check, you should see only the solution given in the paper.