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:
# 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.
The function should then return the modified copy.
Check that your function works properly!
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:
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.
# 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()
# 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$.)
# 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()
# 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.
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:
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.
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?