MAS2008 Scientific Computing: Lab 2
numpy arrays, linear algebra and vectorization


You should complete the tasks below by Wednesday evening of Week 3, and paste your code into the Moodle online test system. You should start working on these tasks during the lab session, and continue afterwards if necessary.

The following notebooks and videos are relevant for this lab:
Numpy arrays View Download
Broadcasting arrays View Download
Linear Algebra View Download
The Mandelbrot set View Download
Universal Functions View Download

Task 1: Row sum norm and condition numbers

You will be asked to upload code for this task to the online test system. The test will allow you to use some numpy functions but not others. When developing your code (but not when submitting it) you should include the following import statement at the top:

from numpy import max, sum, abs, array, zeros
from numpy.linalg import inv
This will allow you to write max(a) for the maximum of the array a. (It would be more usual to have import numpy as np at the top and to write np.max(a), but that will not work in this context.) The same applies for the functions sum, abs, array, zeros and inv. You should not use any other numpy functions in your code for this task, and you should remove the import statement from the code that you submit.

The row sum norm of a matrix $A$ is defined as follows: For each row, we add up the absolute values of the elements in that row, and then take the maximum of these sums. For example, in the $2\times 2$ case, the row-sum norm of a matrix $\begin{bmatrix} a&b\\ c&d\end{bmatrix}$ is $\max(|a|+|b|,|c|+|d|)$.

The condition number of $A$ is the product of the row sum norm of $A$ and the row sum norm of $A^{-1}$. This is important in numerical linear algebra, because calculations with matrices that have a large condition number can be unstable. A well-known example is the Hilbert matrix $H_n$, which has entries $1/(i+j+1)$ for $0\leq i,j\lt n$. For example: \[ H_4 = \begin{bmatrix} 1 & 1/2 & 1/3 & 1/4 \\ 1/2 & 1/3 & 1/4 & 1/5 \\ 1/3 & 1/4 & 1/5 & 1/6 \\ 1/4 & 1/5 & 1/6 & 1/7 \end{bmatrix}. \]

Task 2: Onion matrices

In this task and the following ones, you can use any numpy functions that you like. You should revert to the usual setup: put import numpy as np at the top of your code, and then write np.max(), np.array() and so on.

Given a vector $v=[a,b,c]$ of length $3$ we can construct an "onion matrix" as shown below: \[ \begin{bmatrix} a & a & a & a & a \\ a & b & b & b & a \\ a & b & c & b & a \\ a & b & b & b & a \\ a & a & a & a & a \end{bmatrix} \] This is a $5\times 5$ matrix with a ring of $a$s around the edge, and a ring of $b$s inside that, and a single $c$ in the centre. There is an obvious way to generalise this for vectors of any length. Write a Python function onion(v) which returns the onion matrix corresponding to v, as a numpy array. You can assume that numpy has been imported as np.

You should work out a few examples by hand and check that your code gives the right answer. When you are sure that your code is correct, you should again paste it into the online test on Moodle.

Task 3: Moving average

Fix a number $n\gt 0$. We can define a sequence $a_0,a_1,a_2,\dotsc$ as follows: for $k\lt n$ we have $a_k=k$, and when $k\geq n$ the number $a_k$ is the average of the previous $n$ numbers in the sequence. For example, when $n=3$ the sequence starts like this: \[ a_0=0,\; a_1 = 1,\; a_2=2,\; a_3=(0+1+2)/3 = 1,\; a_4=(1+2+1)/3 = 4/3,\dotsc. \] Write a Python function a(n,k) that calculates $a_k$ for the given value of $n$. Any correct function will be accepted, but ideally you should try to write a function that uses matrix methods and is efficiently vectorized. You can assume that numpy has been imported as np. When you are sure that your code is correct, you should again paste it into the online test on Moodle.

Task 4: Matrix exponentials

Let $A$ be a square matrix, of shape $n\times n$ say. The exponential of $A$ is defined to be $\exp(A)=\sum_{k=0}^\infty A^k/k!$. Here $A^0$ should be interpreted as the identity matrix of size $n$, and $A^1=A$, and in general $A^k$ is the matrix product of $A$ with itself $k$ times. (Just as with the exponential function for real numbers, there are many applications for this.) A reasonable way to calculate the exponential is as follows. We keep track of a number $k$, a term $T$ (which will be $A^k/k!$) and a sum $S$. We start with $k=0$ and $T=I$ and $S=I$. At each step we increase $k$ by $1$, then replace $T$ by $AT/k$, to get the next term, and then add $T$ to $S$. We stop when $T$ is small enough, say when np.max(np.abs(T)) is less than $10^{-9}$.
(This task is not included in the online test.)