Due Date: February 4, 8:59pm

Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html

Please hand in 2 files:

- Python File containing all your code, named
`csc338_a2.py`

. If you are using Jupyter Notebook to complete the work, your notebook can be exported as a .py file (File -> Download As -> Python). Your code will be auto-graded using Python 3.6, so please make sure that your code runs. - PDF file named
`csc338_a2_written.pdf`

containing your solutions to the written parts of the assignment, which includes Q2(e,f), Q3(c,d), Q4(b), Q5(e), and Q6. Your solution can be hand-written, but must be legible. Graders may deduct marks for illegible or poorly presented solutions.

Submit the assignment on **MarkUs** by 9pm on the due date. See the syllabus for the course policy regarding late assignments. All assignments must be done individually.

In [1]:

```
import math
import numpy as np
import matplotlib.pyplot as plt
```

Write a function `solve_diagonal(A, b)`

that takes a diagonal numpy matrix `A`

of size $n \times n$, and a numpy array `b`

of size $n$, and returns the solution `x`

to the problem `Ax=b`

.

Hint: You may find the numpy function `np.diag`

helpful.

In [2]:

```
def solve_diagonal(A, b):
"""Return a vector x with np.matmul(A, x) == b, where
* A is an nxn numpy matrix that is diagonal (A[i,j] == 0 for i != j)
* b is an nx1 numpy vector
>>> A = np.array([[2., 0.],
[0., 3.]])
>>> b = np.array([1., 9.])
>>> solve_diagonal(A, b)
array([0.5, 3.])
"""
```

In [102]:

```
A = np.array([[2., 0., 0.],
[1., 1., 0.],
[0., 1., 2.]])
b = np.array([4., 5., 4.])
q2_x = np.zeros_like(b) # todo
```

In [4]:

```
# Code from tutorial 3
def backward_substitution(A, b):
"""Return a vector x with np.matmul(A, x) == b, where
* A is an nxn numpy matrix that is upper-triangular and nonsingular
* b is an nx1 numpy vector
"""
n = A.shape[0]
x = np.zeros_like(b, dtype=np.float)
for i in range(n-1, -1, -1):
s = 0
for j in range(n-1, i, -1):
s += A[i,j] * x[j]
x[i] = (b[i] - s) / A[i,i]
return x
def eliminate(A, b, k):
"""Eliminate the k-th row of A, in the system np.matmul(A, x) == b,
so that A[i, k] = 0 for i < k. The elimination is done in place."""
n = A.shape[0]
for i in range(k + 1, n):
m = A[i, k] / A[k, k]
for j in range(k, n):
A[i, j] = A[i, j] - m * A[k, j]
b[i] = b[i] - m * b[k]
def gauss_elimination(A, b):
"""Return a vector x with np.matmul(A, x) == b using
the Gauss Elimination algorithm, without partial pivoting."""
for k in range(A.shape[0] - 1):
eliminate(A, b, k)
x = backward_substitution(A, b)
return x
```

For the values of `A`

and `b`

below,
show the updated value of `A`

and `b`

after each elimination step
when solving the linear system $Ax=b$.

Save those values in arrays `q2_A`

and `q2_b`

so that `q2_A[0]`

shows the elimination
of the first column below the diagonal, `q2_A[1]`

shows the elimination
of the second column below the diagonal, and so on. Similarly, `q2_b[i]`

should show the
corresponding value of `b`

the same number of eliminations.

Hint: You can do this question by hand, or write a different version of the `gauss_elimination`

function. You may also find the function `np.copy`

helpful.

In [105]:

```
A = np.array([[2., 0., 0.],
[1., 1., 0.],
[0., 1., 2.]])
b = np.array([4., 5., 4.])
q3_A = []
q3_b = []
```

Write a helper function `partial_pivot`

that performs partial pivoting on `A`

at column `k`

.
Your helper function will be called from the function `gauss_elimination_partial_pivot`

.

In [131]:

```
def partial_pivot(A, b, k):
"""Perform partial pivoting for column k. That is, swap row k
with row j > k so that the new element at A[k,k] is the largest
amongst all other values in column k below the diagonal.
This function should modify A and b in place, and return None.
"""
def gauss_elimination_partial_pivot(A, b):
"""Return a vector x with np.matmul(A, x) == b using
the Gauss Elimination algorithm, with partial pivoting."""
for k in range(A.shape[0] - 1):
partial_pivot(A, b, k)
eliminate(A, b, k)
x = backward_substitution(A, b)
return x
```

Solve the system $Ax = b$ for the values of $A$ and $b$ below using `gauss_elimination`

and then `gauss_elimination_partial_pivot`

.

Save the solution you obtain from `gauss_elimination`

in the variable `q2d_nopivot`

.
Save the solution you obtain from `gauss_elimination_partial_pivot`

in the variable `q2d_pivot`

.

Which is the correct answer? Save the correct answer in the variable `q2d_correct`

.

Note: Make sure your code accounts for the fact that both function `gauss_elimination_partial_pivot`

and `gauss_elimination`

modify their parameters.

In [7]:

```
e = pow(2, -100)
A = np.array([[e, 1],
[1, 1]])
b = np.array([1+e, 2])
q2d_nopivot = None
q2d_pivot = None
q2d_correct = None
```

Why do we obtain different results from the two algorithms? Include your answer in your pdf writeup. Be specific.

Is the matrix $A$ from part (d) well-conditioned or ill-conditioned? Include your answer and justifications in your pdf writeup. Be specific.

Consider the following matrices `M1`

, `M2`

, and `M3`

. Compute each of their $L_1$, $L_2$, and $L_\infty$ norms. Save the results in the variables below.

For the $L_2$ norm you may find the function `np.linalg.norm`

helpful. You can compute the $L_1$ and $L_\infty$ norms either by hand or write a function.

In [100]:

```
M1 = np.array([[3., 0.],
[-4., 2.]])
M2 = np.array([[2., -2., 0.3],
[0.5, 1., 0.9],
[-4., -2., 5]])
M3 = np.array([[0.2, -0.2],
[1.0, 0.2]])
# for matrix M
M1_l_1 = 0
M1_l_2 = 0
M1_l_infty = 0
M2_l_1 = 0
M2_l_2 = 0
M2_l_infty = 0
M3_l_1 = 0
M3_l_2 = 0
M3_l_infty = 0
```

Show that $||A||_1 = max_j \sum_{i=1}^m |a_{ij}|$. Include your solution in your pdf writeup.

Is it true that for a vector, $||v||_1 \le ||v||_\infty$? What about for a matrix: is it true that for a matrix, $||M||_1 \le ||M||_\infty$? Include your solution and justification in your pdf writeup.

Assume $x$ and $y$ are 2-vectors. Is it possible to have $||x||_1 > ||y||_1$ but $||x||_2 < ||y||_2$?

If so, save an example of such vector in the variables `q5_x`

and `q5_y`

. Otherwise, set both
variables to the Python value `None`

.

In [103]:

```
q5_x = np.array([0.0, 0.0])
q5_y = None
```

Classify each of the following matrices `A1`

, `A2`

, `A3`

and `A4`

,
as well-conditioned or ill-conditioned.

You should do this question by hand instead of writing a function to compute condition numbers.

Save the classifications in a Python array called `conditioning`

. Each item of the array
should be either the string "well" or the string "ill".

In [132]:

```
A1 = np.array([[3, 0],
[0, pow(3, -30)]])
A2 = np.array([[pow(3, 20), 0],
[0, pow(5, 50)]])
A3 = np.array([[pow(4, -40), 0],
[0, pow(2, -80)]])
A4 = np.array([[pow(5, -17), pow(5, -16)],
[pow(5, -18), pow(5, -17)]])
# whether A1, A2, A3, A4 is "well" or "ill" conditioned
conditioning = [None, None, None, None]
```

Suppose that $A$ is well-conditioned. Is $A^2$ also well-conditioned? Why or why not? Include your answer and justifications in your pdf writeup. Be specific.

In [133]:

```
def forward_substitution(A, b):
"""Return a vector x with np.matmul(A, x) == b, where
* A is an nxn numpy matrix that is lower-triangular and nonsingular
* b is a nx1 numpy vector
>>> A = np.array([[2., 0.],
[1., -2.]])
>>> b = np.array([1., 2.])
>>> forward_substitution(A, b)
array([ 0.5 , -0.75])
"""
```

Write a function `elementary_elimination_matrix`

that returns the
elements below the $k$-th diagonal in the $k$-th elementary elimination
matrix. These values corresponds to the values of [-$m_{k+1}$, ... - $m_n$]
from your notes.

Since Python indices begin at zero, the first elementary elimination matrix is for $k=0$.

Do not perform any pivoting.

You may assume that `A[i,j] = 0 for i > j, j < k`

In [107]:

```
def elementary_elimination_matrix(A, k):
"""Return the elements below the k-th diagonal of the
k-th elementary elimination matrix.
Precondigion: A[i,j] = 0 for i > j, j < k
>>> A = np.array([[2., 0., 1.],
[1., 1., 0.],
[2., 1., 2.]])
>>> elementary_elimination_matrix(A, 0)
[-0.5, -1.0]
>>> A = np.array([[2., 0., 1.],
[0., 1., 0.],
[0., 1., 2.]])
>>> elementary_elimination_matrix(A, 1)
[-1.0]
"""
```

Write a function `lu_factorize`

that factors a matrix A into its upper and lower
triangular components. Use the function `elementary_elimination_matrix`

as a helper.

In [74]:

```
def lu_factorize(A):
"""Return two matrices L and U, where
* L is lower triangular
* U is upper triangular
* and np.matmul(L, U) == A
>>> A = np.array([[2., 0., 1.],
[1., 1., 0.],
[2., 1., 2.]])
>>> L, U = lu_factorize(A)
>>> L
array([[1. , 0. , 0. ],
[0.5, 1. , 0. ],
[1. , 1. , 1. ]])
>>> U
array([[ 2. , 0. , 1. ],
[ 0. , 1. , -0.5],
[ 0. , 0. , 1.5]])
"""
```

Write a function `solve_lu`

that solves a linear system Ax=b by

```
* factoring $A = LU$ (using the `lu_factorize` function)
* solving $Ly = b$ using forward substitution (using the `forward_substitution` function)
* solving $Ux = y$ using backward substitution (using the `backward_substitution` function)
```

In [130]:

```
def solve_lu(A, b):
"""Return a vector x with np.matmul(A, x) == b using
LU factorization, without partial pivoting.
"""
```

Prove that the matrix

```
A = [[0 1]
[1 0]]
```

has no $LU$ facotirzation. No lower triangular matrix $L$ and upper triangular matrix $U$ such that $A = LU$.

Include your proof in your pdf writeup.

Write a function `invert_matrix`

that takes an $n \times n$ matrix `A`

, and computes its inverse by solving $n$
systems of linear equations of the form $Ax = b$.

You can use any code that you wrote in this assignment. You **will** be graded for efficiency.

Include your code in **both** your `.py`

file **and** your pdf writeup.

In [125]:

```
def invert_matrix(A):
"""Return the inverse of the nxn matrix A by
solving n systems of linear equations of the form
Ax = b.
>>> A = np.array([[0.5, 0.],
[-1., 2.]])
>>> invert_matrix(A)
array([[ 2. , -0. ],
[ 1. , 0.5]])
"""
```