Homework 6

INSERT YOUR NAME HERE

Directions: Add work to this notebook to solve the problems below.

Be careful to name the objects you create as described in the problem. Some problems may include some tests you can run to test your code.

This assignment will not be collected.

Some problems are derived from problems in these books:

  • LL: Programming for Computations - Python by Svein Linge and Hans Petter Langtangen, 2nd edition.
  • L: A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
  • TAK: Applied Scientific Computing With Python by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh.

Standard imports

In [ ]:
import numpy as np

1. Manipulating matrices 1

Write a function manipulate1(M), which takes as input an $m \times n$ matrix with both $m \geq 2$ and $n \geq 2$. The function should make a copy of $M$ (so that the original matrix does not change). Then the function should perform the following operations on the copy in the order below.

  1. It should double the entries in column $1$ (the second column).
  2. It should switch rows $0$ and $1$ (the first and second rows).
  3. It should add three times column $0$ to column $1$.

The function should then return the modified copy.

Check that your function works properly!

In [ ]:

2. Manipulating matrices 2

Write a function manipulate2(m) which takes as input a $2$-dimensional numpy array m representing a matrix $M$. The function should make a copy of m (so that m is not changed). Then the following operations should be performed, in order:

  • All diagonal entries should be doubled.
  • The top right and bottom left entries should be swapped.
  • The entries which are strictly closer to the bottom right corner than any other corner should be negated.

For example, the if the following matrix is input for manipulate2(m) $$\left(\begin{matrix} 2 & 3 & 4 \\ 1 & -1 & 3 \\ 8 & 5 & 11 \end{matrix}\right)$$ then the following matrix should be output: $$\left(\begin{matrix} 4 & 3 & 8 \\ 1 & -2 & 3 \\ 4 & 5 & -22 \end{matrix}\right).$$

In [ ]:

Tests for your code:

In [ ]:
# Test suggested in the problem
m_in = np.array( [ [ 2,  3, 4 ],
                   [ 1, -1, 3 ],
                   [ 8,  5, 11] ] )
m_out = np.array( [ [ 4,  3,  8 ],
                    [ 1, -2,  3 ],
                    [ 4,  5, -22] ] )
assert (manipulate2(m_in) == m_out).all()

3. Manipulating matrices 3

Write a function manipulate3(m) which takes as input a $2$-dimensional numpy array m representing a matrix $M$. The function should return the matrix $$\left(\begin{array}{r|r} 4I & 3 M \\ \hline M^T & M^T M \end{array}\right).$$ Here, $M^T$ represents the transpose of $M$. The matrix $I$ represents the identity matrix whose sizes are chosen to make the above block matrix form well-defined. The matrix should be returned as a $2$-dimensional numpy array.

In [ ]:

4. Inverting a matrix

We have seen that you can solve the matrix equation $A x = b$ using Gaussian elimination, where $A$ is an invertible square matrix. We didn't quite do it this way, but row reduction can turn the block matrix $[A | b]$ into $[I | x]$, where $x$ is the solution.

You can obtain the inverse of a matrix $A$ in this way. Applying row reduction to the block matrix $[A | I]$ leads to the block matrix $[I | A^{-1}]$. Write a function inverse(A) which takes as input a invertible square matrix of floats, $A$, and returns the inverse $A^{-1}$. You must find this inverse using row reduction. Recall that allowed moves in row reduction are:

  • Swapping two rows.
  • Adding a multiple of one row to another row.
  • Multiplying a row by a scalar.

Remarks: You need to use pivoting to deal with the fact that some pivots might be zero. You can check your work easily. If $M$ is an invertible $n \times n$ matrix, then the following should return a small matrix:

np.abs( M @ inverse(M) - np.identity(n) )

Some tests are provided below.

In [ ]:

5. Tridiagonal matrices

Write a function tridiagonal_matrix(a,b,c) which takes as input a vector $a$ with some length $n$ and vectors $b$ and $c$ with length $n-1$. The function should return the matrix $$\left[\begin{array}{rrrrrr} a_0 & b_0 & & & & \\ c_0 & a_1 & b_1 & & & \\ & c_1 & a_2 & b_2 & & \\ & & \ddots & \ddots & \ddots & \\ & & & c_{n-3} & a_{n-2} & b_{n-2} \\ & & & & c_{n-2} & a_{n-1} \end{array}\right].$$

(This matrix appears in equation (4.10) of TAK on page 94, but I have shifted the indices to match Python conventions that we start indexing at zero rather than $1$.)

In [ ]:

Tests for your code:

In [ ]:
assert \
    ( tridiagonal_matrix(np.arange(3),np.arange(3,5),np.arange(5,7)) == \
      np.array([[0., 3., 0.],
                [5., 1., 4.],
                [0., 6., 2.]]) ).all(), \
    "Error. The matrix returned is incorrect."

6. Solving tridiagonal systems

Read about tridiagonal systems in TAK § 4.2.2.

Write a function tridiagonal_solve(a, b, c, r) to solve a tridiagonal system using Gauss elimination. The program should take as input four vectors: A vector $a$ in ${\mathbb R}^n$ for some $n$, vectors $b$ and $c$ in ${\mathbb R}^{n-1}$ and $r \in {\mathbb R}^n$. The function should return a vector $x$ that solves the equaion $Mx=r$, where $M$ is the $n \times n$ matrix with $a$ on the diagonal, $b$ above the diagonal, and $c$ below the diagonal. See equation (4.10) on page 94 of TAK.

Test that your function works on a $50 × 50$ system with any diagonal entries of your choice and a right hand side so that the exact solution vector is the vector of all ones.

Hint: The algorithm is spelled out as algorithm 6 on page 94 of TAK, but be wary of the indexing changed mentioned in the previous problem. You may assume that there are no zero pivots ($a_{i} \neq 0$) as you row reduce. This is what TAK assumes. You can check your work using the prior problem. For reasonable inputs $a$, $b$, $c$ and $r$, the following test should pass:

M = tridiagonal_matrix(a, b, c)
x = tridiagonal_solve(a, b, c, r)
assert (np.abs(M @ x - r) < 10**-8).all()
In [ ]:

Tests for your code:

In [ ]:
a = np.array([2., 2., 2.])
b = np.array([1., 1.])
c = np.array([1., 1.])
r = np.array([4., 8., 8.])
M = tridiagonal_matrix(a, b, c)
x = tridiagonal_solve(a, b, c, r)
assert (np.abs(M @ x - r) < 10**-8).all()
print("Passed!")
# The solution should be [1 2 3]
print(x)