Student's Name: PLEASE INSERT YOUR NAME HERE

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

Check your work. I give some tests that should be passed if you get it correct.

**Your notebook should run without printing errors and without user input. If this is not the case, points will be deducted from your grade.**

Problem Sources:

**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.

In [ ]:

```
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv
```

Write a function `manipulate(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.

- It should double the entries in column $1$ (the second column).
- It should switch rows $0$ and $1$ (the first and second rows).
- 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 [ ]:

```
```

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 [ ]:

```
```

In [ ]:

```
# Test on the matrix A
A = np.array([[2.,1.],
[1.,1.]])
# The actual inverse:
Ai = np.array([[ 1., -1.],
[-1., 2.]])
# Test if the returned matrix is close to the actual inverse
assert(inverse(A)-Ai < 10**-8).all()
```

In [ ]:

```
# Test on some random 2x2 matrices with entries in [-1,1]
num_tests = 3
for i in range(num_tests):
A = 2*np.random.random(4).reshape(2,2) - 1
print("Testing on \n{}".format(A))
Ai = inverse(A)
assert (np.abs(A @ Ai - np.identity(2)) < 10**-8).all(), "Test failed."
print("Passed!\n")
```

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 [ ]:

```
```

In [ ]:

```
# Tests:
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."
```

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 [ ]:

```
```

In [ ]:

```
# One test
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)
```

(Before starting this problem read Example 7, which starts at the bottom of page 94 of TAK. You may want to go back and read Example 1 on page 36.)

Consider solving the Bessel equation $$x^2 y'' + x y' + (x^2 - 1) y = 0$$ subject to the boundary conditions $y(1)=1$ and $y(15)=0$ using numerical derivatives with $N = 280$ steps between $[1, 15]$. You should begin by rewriting the differential equation in the form $$y'' + \frac{1}{x} y' + (1-\frac{1}{x^2}) y = 0.$$ Next, use finite difference approximations to the derivatives and collect like terms to obtain a tridiagonal linear system for the unknown discrete values of y.

Use the function `tridiagonal_solve`

that you wrote to find the approximate solution to the system. Graph the solution.

(This is problem 11 from TAK § 4.2.)

*First hint:* 280 points dividing $[1, 15]$ into equal sized intervals are computed below.

In [ ]:

```
N = 280
x = np.linspace(1,15,N+2)[1:-1]
assert len(x)==N
print("x[0] = {}".format(x[0]))
print("x[N-1] = {}".format(x[N-1]))
```

Our step size $h$ is therefore:

In [ ]:

```
h = x[0] - 1
print(h)
```

We are looking to find a numpy array $y$ which can be obtained by multiplying by a tridiagonal matrix to get some vector we know. Then we can use `tridiagonal_solve`

to find $y$. Here the entries of $y$, $y_i$, should be given by $y_i=y(x_i)$, where $x_i$ is the numpy array above. Thus typically $y(x_i+h)=y_{i+1}$ and $y(x_i-h)=y_{i-1}$ (but there are special cases at the end of the arrays).

Observe that $$\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] \left(\begin{array}{r} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-2} \\ y_{N-1}\end{array}\right)= \left(\begin{array}{r} a_0 y_0 + b_0 y_1 \\ c_0 y_0 + a_1 y_1 + b_1 y_2 \\ c_1 y_1 + a_2 y_2 + b_2 y_3 \\ \vdots \\ c_{N-3} y_{N-3} + a_{N-2} y_{N-2} + b_{N-2} y_{N-1} \\ c_{N-2} y_{N-2} + a_{N-1} y_{N-1} \end{array}\right)=r. $$ You can find the values of $a$, $b$, $c$ and $r$ by analyzing the differential equation. You want to replace $y'$ by the two-point approximation of the dertivative (with $h$ as computed above), and $y''$ with the three point approximation of the derivative.

In [ ]:

```
```

Solve $y'' + 2x y + 2y = 3e^{−x} − 2xe^{−x}$ subject to $y(−1) = e + 3/e$, $y (1) = 4/e$ using the second-order finite difference method as in the above problem. Consider $N = 10, 20, 40$ in $[−1, 1]$. Plot these solutions and the true solution $y = e^{−x} + 3e^{−x^2}$. Also, on a separate set of axes, plot the three error curves. Estimate the order of the truncation error for this method. Is it what you would expect?

In [ ]:

```
```