<jupyter_start><jupyter_text>Turn in an image (e.g., screenshot) or PDF copy of any code that is part of your answer. Make sure all images and PDF pages are properly rotated. Make sure that all pages are clearly visible. Tips: Use the document scanner function on your smart phone to take better page "scans" using your camera. Make sure your screen is not shifted toward warmer colours (some devices filter blue light at night) giving it a dim and orange appearance.   Q1 ADerive a method for computing the determinant of a matrix $A\in \mathbb{R}^{n\times n}$ using Gaussian elimination with partial pivoting. **Hint: use the fact that a permutation matrix is an orthogonal matrix.**------------------------------------------------------------------------------For Gaussian elimination with partial pivoting we have $$ P A = LU, $$where $P$ is a permutation matrix. Since the permutation matrix is an orthogonal matrix, it is also nonsingular. Taking the determinant to both sides of the above equation yields $\det(P) \det(A) = \det(L)\det(U)$. Since $\det(L) = 1$ we have that $\det(A) = \det(U)/\det(P)$. If we know $\det(P)$ then our formula is given by $\det(A) = \frac{1}{\det(P)} \prod_{i=1}^n u_{ii}$. This last follows from the fact that $U$ is upper triangular.Since $P$ is an orthogonal matrix we have that $I = P^T P$, and it follows that $1 = \det(I) = \det(P^T P) = \det(P^T)\det(P) = \det(P)^2$. Hence $\det(P) = \pm 1$. to determine the sign of the determinant, we define $m$ to be the number of row swaps used in Gauassian elimination with partial pivoting, and we define $\sigma = (-1)^m$. The matrix $P$ is rarely formed explicitly. Instead a vector of row permutations is stored and one can keep track of the number of row swaps that are used. One can show that $\det(P) = \sigma$. This last step requires no additional floating point operations to compute, since the information is obtained from the decomposition without extra work.Our final formula becomes$$\det(A) = \sigma\prod_{i=1}^n u_{ii}.$$ BAssume that the LU decomposition has already been computed. Show that the method for computing the determinant requires $n$ floating point operations.------------------------------------------------------------------------------The determinant formula in part A requires exactly $n$ multiplications.  Q2Let $b + \delta b$ be a perturbation of the vector $b\neq 0$ and let $x$ and $\delta x$ be such that $Ax = b$ and $A(x + \delta x) = b + \delta b$, where $A$ is a given nonsingular matrix. Show that $$\frac{\Vert \delta x \Vert }{\Vert x \Vert } \leq \kappa(A) \frac{\Vert \delta b \Vert }{\Vert b \Vert }.$$------------------------------------------------------------------------------From $A(x + \delta x) = b + \delta b$ we use $Ax = b$ and solve for $\delta x$ to get $ A\delta x = \delta b $.Since $A$ is invertible, we have $\delta x = A^{-1}\delta b$. Taking norms to both sides yields $\Vert \delta x \Vert =  \Vert A^{-1}\delta b \Vert$. Using the submultiplicative property, we get $\Vert \delta x \Vert \leq  \Vert A^{-1}\Vert  \Vert \delta b \Vert$. Since $A$ is nonsingular and $b \neq 0$, we have that $\Vert x \Vert \neq 0$. We can multiply and divide by $\Vert x \Vert$ to get $\frac{\Vert \delta x \Vert}{\Vert x \Vert} \leq  \Vert A^{-1}\Vert  \frac{\Vert \delta b \Vert}{\Vert x \Vert}$. Using $\Vert b \Vert \leq \Vert A \Vert \Vert x \Vert$ and $\Vert A \Vert \neq 0$ (since $A$ is nonsingular) we have $\frac{\Vert \delta x \Vert}{\Vert x \Vert} \leq  \Vert A^{-1}\Vert  \frac{\Vert \delta b \Vert}{\Vert b \Vert / \Vert A \Vert }$. The result follows after substitution of $\kappa(A) = \Vert A \Vert \Vert A^{-1} \Vert$.  Q3The book states without proof that if $\Vert \delta A \Vert < \frac{1}{\Vert A^{-1}\Vert}$ then the bound on the relative error can be written as$$  \frac{\Vert x - \hat{x} \Vert}{\Vert x \Vert} \leq \frac{\kappa(A)}{1 - \kappa(A)\frac{\Vert \delta A\Vert}{\Vert A\Vert}}  \left(\frac{\Vert \delta b \Vert}{\Vert b\Vert} + \frac{\Vert\delta A \Vert}{\Vert A\Vert} \right). $$In what follows, let $A$, $B$, and $\delta A$ be real $n\times n$ matrices. We saw an incomplete proof in lecture that assumed $A + \delta A$ is nonsingular. In the following, you will show that if $\Vert \delta A \Vert < \frac{1}{\Vert A^{-1}\Vert}$ then $A + \delta A$ is nonsingular. AShow that if the spectral radius $\rho(A) < 1$ then the matrix $A - I$ is nonsingular. **Hint: use the definition of the spectral radius and the characteristic equation, $\det(A-\lambda I) = 0$, for the eigenvalues of $A$.**---------------------------------------------------------------------------------If $\rho(A) < 1$ then by definition all of the eigenvalues $\vert \lambda \vert < 1$, which means that there are no eigenvalues on the unit circle in the complex plane. It follows that $\det(A \pm I) \neq 0$. Hence $A - I$ is nonsingular. BShow that if $A$ is nonsingular and $\Vert A - B\Vert  < \frac{1}{\Vert A^{-1} \Vert}$ then $B$ is nonsingular. **Hint: use $B = A[I - A^{-1}(A - B)]$ and part A**---------------------------------------------------------------------------------We can write $B = A[I - A^{-1}(A - B)]$. It follows that since $A$ is nonsingular then $B$ is nonsingular if $I- A^{-1}(A - B)$ is nonsingular. From part A, we know that $I- A^{-1}(A - B)$ is nonsingular if $\rho( A^{-1}(A - B)) < 1$. To show this we use the fact that for any square matrix $M$ that $\rho(M) \leq \Vert M\Vert$, and we use the assumption $\Vert A - B\Vert  < \frac{1}{\Vert A^{-1} \Vert}$. We have that$$\rho(A^{-1}(A - B)) \leq \Vert A^{-1}(A - B) \Vert \leq \Vert A^{-1}\Vert \Vert (A - B) \Vert  < 1.$$ CShow that if $A$ is nonsingular and $\Vert \delta A\Vert  < \frac{1}{\Vert A^{-1} \Vert}$, then $A + \delta A$ is nonsingular.---------------------------------------------------------------------------------Let $B = A + \delta A$ and apply part B.  Q4Tridiagonal matrices appear often. Due to their simple structure, it is possible to significantly speed up the computation of the LU decomposition. AImplement Gaussian elimination for computing $A = LU$ of the form\begin{equation}\begin{bmatrix}a_1 & c_1 & & & & \\b_2 & a_2 & c_2& & & \\ & b_3& a_3& c_3& & \\ & &\ddots & \ddots & \ddots & \\ & & & \ddots & \ddots& c_{n-1}\\ & & & &b_n & a_n\end{bmatrix} =\begin{bmatrix}1 & & & & & \\l_2 & 1& & & & \\ & l_3 & 1& & & \\ & & \ddots & \ddots & & \\ & & & \ddots& \ddots& \\ & & & & l_n & 1\end{bmatrix}\begin{bmatrix} u_1& c_1& & & & \\ & u_2& c_2 & & & \\ & & u_3& c_3 & & \\ & & & \ddots& \ddots& \\ & & & & \ddots& c_{n-1}\\ & & & & & u_n\end{bmatrix}\end{equation}Your function should take three input arguments: the vectors `a`, `b`, and `c` containing the diagonals of the matrix $A$. It should return two vectors `l` and `u` containing the elements from the above $LU$ decomposition (ie the vectors $l$ and $u$ have entries $l_i$ and $u_i$).  SolutionWe usually index vectors and matrices in linear algebra starting with the value 1. For example, the first value of a vector $x$ is $x_1$. In Python and most programming languages, indexing of array starts with the value 0. For example, the first value stored in an array `x` is `x[0]`. For this reason, our formulas and our code will use different indices. Let $i$ be the index of the vector $x$, and let `k` be the index for the array `x`. We set $i = $ `k+1`. Then $x_i = $ `x[i-1]` and `x[k]` $=x_{k+1}$.  Version 1In this version, I pad the arrays `b` and `l` with a leading `nan` so that the indices line up. This way, $b_2$ corresponds to `b[1]` and $l_2$ corresponds to `l[1]`. Note that this is consistent with the other vectors; i.e., $a_1$ corresponds to `a[0]`.<jupyter_code>def tridiag_solve_version1(a, b, c):
    n = a.size
    u = zeros(n)
    u[0] = a[0]
    l = zeros(n)
    l[0] = nan
    for i in arange(2, n+1):
        k = i-1
        l[k] = b[k]/u[k-1]
        u[k] = a[k] - l[k]*c[k-1]
    return l, u<jupyter_output><empty_output><jupyter_text>Version 2In this version, I do not pad the arrays `b` and `l` with a leading entry to align the indices. Instead, I adjust the formulas so that $u_2$ and $l_2$ correspond to `u[0]` and `l[0]`, respectively. In other words, there is a different relationship between $i$ and the index for `u` and `l`. Let `p` $=i-2$ index these vectors, but not the others (for the others we have `k` $=i-1$).<jupyter_code>def tridiag_solve_version2(a, b, c):
    n = a.size
    u = zeros(n)
    u[0] = a[0]
    l = zeros(n-1)
    for i in arange(2, n+1):
        k = i-1
        p = i-2
        l[p] = b[p]/u[k-1]
        u[k] = a[k] - l[p]*c[k-1]
    return l, u<jupyter_output><empty_output><jupyter_text>BUse your function from part A to compute the LU decomposition of $$ A = \begin{bmatrix}1 & -\frac{1}{2} & & & & \\-\frac{2}{2} & 2& -\frac{2}{2}& & & \\ & -\frac{3}{2}& 3& -\frac{3}{2}& & \\ & & \ddots & \ddots& \ddots& \\ & & & \ddots& \ddots& -\frac{10-1}{2} \\ & & & & -\frac{10}{2}& 10\end{bmatrix}.$$Include a print out of the elements of `l` and `u` rounded to three decimal places.  Solution Version 1In this version, I pad the arrays `b` and `l` with a leading `nan` so that the indices line up. This way, $b_2$ corresponds to `b[1]` and $l_2$ corresponds to `l[1]`. Note that this is consistent with the other vectors; i.e., $a_1$ corresponds to `a[0]`.<jupyter_code>n = 10
a = arange(n) + 1.
b = -(arange(n)-1 + 2.)/2.
b[0] = nan
c = -(arange(n-1) + 1.)/2.


l, u = tridiag_solve_version1(a, b, c)

print(around(l, 3))
print(around(u, 3))
## check 
a_check = c*l[1:] + u[1:]
## output should be 2, 3, 4, ..., 10
print(around(a_check, 3))

b_check = l[1:]*u[:-1]
## output should be -1, -1.5, -2.5, ..., -5
print(around(b_check, 3))
## check 2
from scipy.linalg import lu as lu
A = diag(b[1:], -1) + diag(a) + diag(c, 1)
_, L, U = lu(A)
print(around(L, 2))
print(around(U, 2))
print(norm(L@U - A, 2))<jupyter_output>[[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [-1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. -1.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0. -1.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0. -1.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0. -1.  1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0. -1.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0. -1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0. -1.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0. -1.  1.]]
[[ 1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0. ]
 [ 0.   1.5 -1.   0.   0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   2.  -1.5  0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   0.   2.5 -2.   0.   0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   3.  -2.5  0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.   3.5 -3.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.   0.   4.  -3.5  0.   0. ]
 [ 0.   0.   0.   0.   0.   0.   0.   4.5 -4.   0. ]
 [ 0.   0.   0.   0.   0.   0.   0.   0.   5.  -4.5]
 [ 0.   0.   0.   0.   0.   0.   0.   0.   0.   5.5]]
0.0<jupyter_text>Version 2In this version, I do not pad the arrays `b` and `l` with a leading entry to align the indices. Instead, I adjust the formulas so that $u_2$ and $l_2$ correspond to `u[0]` and `l[0]`, respectively. In other words, there is a different relationship between $i$ and the index for `u` and `l`. Let `p` $=i-2$ index these vectors, but not the others (for the others we have `k` $=i-1$).<jupyter_code>n = 10
a = arange(n) + 1.
b = -(arange(n-1) + 2.)/2.
c = -(arange(n-1) + 1.)/2.


l, u = tridiag_solve_version2(a, b, c)

print(around(l, 3))
print(around(u, 3))
## check 
a_check = c*l + u[1:]
## output should be 2, 3, 4, ..., 10
print(around(a_check, 3))

b_check = l*u[:-1]
## output should be -1, -1.5, -2.5, ..., -5
print(around(b_check, 3))
## check 2
from scipy.linalg import lu as lu
A = diag(b, -1) + diag(a) + diag(c, 1)
_, L, U = lu(A)
print(around(L, 2))
print(around(U, 2))
print(norm(L@U - A, 2))<jupyter_output>[[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [-1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. -1.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0. -1.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0. -1.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0. -1.  1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0. -1.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0. -1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0. -1.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0. -1.  1.]]
[[ 1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0. ]
 [ 0.   1.5 -1.   0.   0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   2.  -1.5  0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   0.   2.5 -2.   0.   0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   3.  -2.5  0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.   3.5 -3.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.   0.   4.  -3.5  0.   0. ]
 [ 0.   0.   0.   0.   0.   0.   0.   4.5 -4.   0. ]
 [ 0.   0.   0.   0.   0.   0.   0.   0.   5.  -4.5]
 [ 0.   0.   0.   0.   0.   0.   0.   0.   0.   5.5]]
0.0
{'fraction_non_alphanumeric': 0.21317252721434723, 'fraction_numerical': 0.05012138773592294, 'mean_word_length': 1.9676969556123636, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '516072', 'license': 'Apache-2.0', 'n_tokens_mistral': 5135, 'n_tokens_neox': 4938, 'n_words': 1931, 'path': 'Homework 7 Solutions.ipynb', 'repo_name': 'newby-jay/MATH381-Fall2021-JupyterNotebooks', 'size': '12769'}