{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 6\n",
"\n",
"Student's Name: PLEASE INSERT YOUR NAME HERE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Directions:__ Add work to this notebook to solve the problems below. \n",
"\n",
"Check your work. I give some tests that should be passed if you get it correct.\n",
"\n",
"**Your notebook should run without printing errors and without user input. If this is not the case, points will be deducted from your grade.**\n",
"\n",
"Problem Sources:\n",
"* __LL__: *Programming for Computations - Python* by Svein Linge and Hans Petter Langtangen, 2nd edition.\n",
"* __L__: *A Primer on Scientific Programming with Python* by Hans Petter Langtangen, 2nd edition.\n",
"* __TAK__: *Applied Scientific Computing With Python* by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Standard imports\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import math as m\n",
"from mpmath import mp, iv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Manipulating matrices\n",
"\n",
"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.\n",
"1. It should double the entries in column $1$ (the second column).\n",
"2. It should switch rows $0$ and $1$ (the first and second rows).\n",
"3. It should add three times column $0$ to column $1$.\n",
"\n",
"The function should then return the modified copy.\n",
"\n",
"*Check that your function works properly!*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Inverting a matrix\n",
"\n",
"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.\n",
"\n",
"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:\n",
"* Swapping two rows. \n",
"* Adding a multiple of one row to another row.\n",
"* Multiplying a row by a scalar.\n",
"\n",
"*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:\n",
"```python\n",
"np.abs( M @ inverse(M) - np.identity(n) )\n",
"```\n",
"Some tests are provided below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test on the matrix A\n",
"A = np.array([[2.,1.],\n",
" [1.,1.]])\n",
"# The actual inverse:\n",
"Ai = np.array([[ 1., -1.],\n",
" [-1., 2.]])\n",
"# Test if the returned matrix is close to the actual inverse\n",
"assert(inverse(A)-Ai < 10**-8).all()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# Test on some random 2x2 matrices with entries in [-1,1]\n",
"num_tests = 3\n",
"for i in range(num_tests):\n",
" A = 2*np.random.random(4).reshape(2,2) - 1\n",
" print(\"Testing on \\n{}\".format(A))\n",
" Ai = inverse(A)\n",
" assert (np.abs(A @ Ai - np.identity(2)) < 10**-8).all(), \"Test failed.\"\n",
" print(\"Passed!\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Tridiagonal matrices\n",
"\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\n",
"$$\\left[\\begin{array}{rrrrrr}\n",
"a_0 & b_0 & & & & \\\\\n",
"c_0 & a_1 & b_1 & & & \\\\\n",
" & c_1 & a_2 & b_2 & & \\\\\n",
" & & \\ddots & \\ddots & \\ddots & \\\\\n",
" & & & c_{n-3} & a_{n-2} & b_{n-2} \\\\\n",
" & & & & c_{n-2} & a_{n-1}\n",
"\\end{array}\\right].$$\n",
"\n",
"(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$.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tests:\n",
"assert \\\n",
" ( tridiagonal_matrix(np.arange(3),np.arange(3,5),np.arange(5,7)) == \\\n",
" np.array([[0., 3., 0.],\n",
" [5., 1., 4.],\n",
" [0., 6., 2.]]) ).all(), \\\n",
" \"Error. The matrix returned is incorrect.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Solving tridiagonal systems\n",
"\n",
"Read about tridiagonal systems in TAK § 4.2.2. \n",
"\n",
"Write a function `tridiagonal_solve(a, b, c, r)` to solve a tridiagonal system using Gauss elimination. \n",
"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. \n",
"\n",
"Test that your function works on a $50 × 50$ system with any diagonal entries of your choice and a right hand\n",
"side so that the exact solution vector is the vector of all ones.\n",
"\n",
"*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:\n",
"```python\n",
"M = tridiagonal_matrix(a, b, c)\n",
"x = tridiagonal_solve(a, b, c, r)\n",
"assert (np.abs(M @ x - r) < 10**-8).all()\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# One test\n",
"a = np.array([2., 2., 2.])\n",
"b = np.array([1., 1.])\n",
"c = np.array([1., 1.])\n",
"r = np.array([4., 8., 8.])\n",
"M = tridiagonal_matrix(a, b, c)\n",
"x = tridiagonal_solve(a, b, c, r)\n",
"assert (np.abs(M @ x - r) < 10**-8).all()\n",
"print(\"Passed!\")\n",
"# The solution should be [1 2 3]\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. Solving a differential equation\n",
"\n",
"(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.)\n",
"\n",
"Consider solving the Bessel equation\n",
"$$x^2 y'' + x y' + (x^2 - 1) y = 0$$ \n",
"subject to the boundary conditions $y(1)=1$ and $y(15)=0$ using numerical derivatives with $N = 280$ steps between $[1, 15]$.\n",
"You should begin by rewriting the differential equation in the form\n",
"$$y'' + \\frac{1}{x} y' + (1-\\frac{1}{x^2}) y = 0.$$\n",
"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. \n",
"\n",
"Use the function `tridiagonal_solve` that you wrote to find the approximate solution to the system. Graph the solution.\n",
"\n",
"(This is problem 11 from TAK § 4.2.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*First hint:* 280 points dividing $[1, 15]$ into equal sized intervals are computed below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"N = 280\n",
"x = np.linspace(1,15,N+2)[1:-1]\n",
"assert len(x)==N\n",
"print(\"x[0] = {}\".format(x[0]))\n",
"print(\"x[N-1] = {}\".format(x[N-1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our step size $h$ is therefore:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"h = x[0] - 1\n",
"print(h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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).\n",
"\n",
"Observe that \n",
"$$\\left[\\begin{array}{rrrrrr}\n",
"a_0 & b_0 & & & & \\\\\n",
"c_0 & a_1 & b_1 & & & \\\\\n",
" & c_1 & a_2 & b_2 & & \\\\\n",
" & & \\ddots & \\ddots & \\ddots & \\\\\n",
" & & & c_{N-3} & a_{N-2} & b_{N-2} \\\\\n",
" & & & & c_{N-2} & a_{N-1}\n",
"\\end{array}\\right]\n",
"\\left(\\begin{array}{r}\n",
"y_0 \\\\\n",
"y_1 \\\\\n",
"y_2 \\\\\n",
"\\vdots \\\\\n",
"y_{N-2} \\\\\n",
"y_{N-1}\\end{array}\\right)=\n",
"\\left(\\begin{array}{r}\n",
"a_0 y_0 + b_0 y_1 \\\\\n",
"c_0 y_0 + a_1 y_1 + b_1 y_2 \\\\\n",
"c_1 y_1 + a_2 y_2 + b_2 y_3 \\\\\n",
"\\vdots \\\\\n",
"c_{N-3} y_{N-3} + a_{N-2} y_{N-2} + b_{N-2} y_{N-1} \\\\\n",
"c_{N-2} y_{N-2} + a_{N-1} y_{N-1}\n",
"\\end{array}\\right)=r.\n",
"$$\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6. Another differential equation\n",
"\n",
"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.\n",
"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.\n",
"Estimate the order of the truncation error for this method. Is it what you would\n",
"expect?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}