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

Instructions

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}. \]

Write a function hilbert_matrix(n) which returns the Hilbert matrix $H_n$.
Wikipedia says that the entries in the Hilbert matrix are $1/(i+j-1)$, which is not the same as our formula above. Is this an error?
Write a function row_sum_norm(A) which returns the row sum norm of a matrix $A$.
Write a function condition_number(A) which returns the condition number of a matrix $A$.
Use web search or an AI assistant to find out how to calculate the condition number using a single numpy function. The relevant function takes two arguments, one of which is the matrix $A$. You may need to search some more or read the documentation to find out what the other argument needs to be.
After adjusting your import statement as necessary, check that your function gives the same answer as the numpy function for some small matrices of your choice, and also for Hilbert matrices of sizes $2$ to $10$.
Paste your functions row_sum_norm(A) and condition_number(A) into the online test on Moodle. Do not include the import statement in the code that you submit. Also, do not include any print statements or test cases; just submit the definitions of the two functions. Your functions should be named exactly as given here, and should take exactly the arguments given here.

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}$.

Write a function matrix_exp(A) which implements this method. You can assume that numpy has been imported as np.
Explain why np.exp(A) is not a solution to the above question. Explain why np.exp(A) gives an answer for non-square matrices, but matrix_exp(A) just gives an error.