{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Laboratory 3: Linear Algebra (Sept. 12, 2017)\n", "\n", "Grace Yung" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Quiz on Matrices\n", "Which matrix matches this system of equations?\n", "
\n", "$$\\begin{array}{lcr}\n", "2x + 3y + 6z &=& 19\\\\\n", "3x + 6y + 9z &=& 21\\\\\n", "x + 5y + 10z &=& 0\n", "\\end{array}$$\n", "
\n", "(A)\n", "$$\\left[ \\begin{array}{ccc}\n", "2 & 3 & 1 \\\\\n", "3 & 6 & 5 \\\\\n", "6 & 9 & 10\\\\\n", "19 & 21 & 0\n", "\\end{array}\n", "\\right]$$\n", "(B)\n", "$$\\left[ \\begin{array}{ccc}\n", "2 & 3 & 6\\\\\n", "3 & 6 & 9\\\\\n", "1 & 5 & 10\n", "\\end{array}\n", "\\right]$$\n", "(C)\n", "$$\\left[ \\begin{array}{ccc|c}\n", "1 & 5 & 10 & 0\\\\\n", "2 & 3 & 6 & 19\\\\\n", "3 & 6 & 9 & 21\n", "\\end{array}\n", "\\right]$$\n", "(D)\n", "$$\\left[ \\begin{array}{ccc|c}\n", "2 & 3 & 6 & -19\\\\\n", "3 & 6 & 9 & -21\\\\\n", "1 & 5 & 10 & 0 \n", "\\end{array}\n", "\\right]$$\n", "\n", "In the following, replace 'x' by 'A', 'B', 'C', or 'D' and run the cell." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\alpha = \\beta$\n", "$$\\int \\alpha dx = \\beta$$\n", "$\\require{ams}$" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "$~~~$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## List of Problems\n", "\n", "- [Problem One](#Problem-One): Pollution Box Model\n", "- [Problem Two](#Problem-Two): Condition number for Dirichlet problem\n", "- [Problem Three](#Problem-Three): Condition number for Neumann problem\n", "- [Problem Four](#Problem-Four): Condition number for periodic problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Objectives\n", "\n", "The object of this lab is to familiarize you with some of the common\n", "techniques used in linear algebra. You can use the Python software\n", "package to solve some of the problems presented in the lab. There are\n", "examples of using the Python commands as you go along in the lab. In\n", "particular, after finishing the lab, you will be able to\n", "\n", "- Define: condition number, ill-conditioned matrix, singular matrix,\n", " LU decomposition, Dirichlet, Neumann and periodic boundary\n", " conditions.\n", "\n", "- Find by hand or using Python: eigenvalues, eigenvectors, transpose,\n", " inverse of a matrix, determinant.\n", "\n", "- Find using Python: condition numbers.\n", "\n", "- Explain: pivoting.\n", "\n", "There is a description of using Numpy and Python for Matrices at the end of the lab. It includes a brief description of how to use the built-in functions introduced in\n", "this lab. Just look for the paw prints:\n", "\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n", "\n", "when you are not sure what\n", "functions to use, and this will lead you to the mini-manual." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "You should have had an introductory course in linear algebra." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "******************************\n", "context imported. Front of path:\n", "/Users/phil/repos/numeric\n", "back of path: /Users/phil/.ipython\n", "******************************\n", "\n", "through /Users/phil/repos/numeric/notebooks/lab3/context.py\n" ] } ], "source": [ "# import the context to find files\n", "import context\n", "# import the quiz script\n", "from numlabs.lab3 import quiz3\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear Systems\n", "\n", "In this section, the concept of a matrix will be reviewed. The basic\n", "operations and methods in solving a linear system are introduced as\n", "well.\n", "\n", "Note: **Whenever you see a term you are not familiar with, you can find\n", "a definition in the [glossary](#Glossary).**\n", "\n", "### What is a Matrix?\n", "\n", "Before going any further on how to solve a linear system, you need to\n", "know what a linear system is. A set of m linear equations with n\n", "unknowns: \n", "\n", "
\n", "(System of Equations)\n", "$$\\begin{array}{ccccccc}\n", "a_{11}x_{1} & + & \\ldots & + & a_{1n}x_{n} & = & b_{1} \\\\\n", "a_{21}x_{1} & + & \\ldots & + & a_{2n}x_{n} & = & b_{2} \\\\\n", " & & \\vdots & & & & \\vdots \\\\\n", "a_{m1}x_{1} & + & \\ldots & + & a_{mn}x_{n} & = & b_{m}\n", "\\end{array}\n", "$$\n", "
\n", "\n", "can be represented as an augmented matrix:\n", "\n", "$$\\left[\n", "\\begin{array}{cc}\n", " \\begin{array}{ccccc}\n", " a_{11} & & \\ldots & & a_{1n} \\\\\n", " \\vdots & & \\ddots & & \\vdots \\\\\n", " a_{m1} & & \\ldots & & a_{mn}\n", " \\end{array}\n", "&\n", " \\left|\n", " \\begin{array}{rc}\n", " & b_{1} \\\\ & \\vdots \\\\ & b_{m}\n", " \\end{array}\n", " \\right.\n", "\\end{array}\n", "\\right]$$\n", "\n", "Column 1 through n of this matrix contain the coefficients $a_{ij}$ of\n", "the unknowns in the set of linear equations. The right most column is\n", "the *augmented column*, which is made up of the coefficients of the\n", "right hand side $b_i$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Quiz on Matrices\n", "Which matrix matches this system of equations?\n", "
\n", "$$\\begin{array}{lcr}\n", "2x + 3y + 6z &=& 19\\\\\n", "3x + 6y + 9z &=& 21\\\\\n", "x + 5y + 10z &=& 0\n", "\\end{array}$$\n", "
\n", "(A)\n", "$$\\left[ \\begin{array}{ccc}\n", "2 & 3 & 1 \\\\\n", "3 & 6 & 5 \\\\\n", "6 & 9 & 10\\\\\n", "19 & 21 & 0\n", "\\end{array}\n", "\\right]$$\n", "(B)\n", "$$\\left[ \\begin{array}{ccc}\n", "2 & 3 & 6\\\\\n", "3 & 6 & 9\\\\\n", "1 & 5 & 10\n", "\\end{array}\n", "\\right]$$\n", "(C)\n", "$$\\left[ \\begin{array}{ccc|c}\n", "1 & 5 & 10 & 0\\\\\n", "2 & 3 & 6 & 19\\\\\n", "3 & 6 & 9 & 21\n", "\\end{array}\n", "\\right]$$\n", "(D)\n", "$$\\left[ \\begin{array}{ccc|c}\n", "2 & 3 & 6 & -19\\\\\n", "3 & 6 & 9 & -21\\\\\n", "1 & 5 & 10 & 0 \n", "\\end{array}\n", "\\right]$$\n", "\n", "In the following, replace 'x' by 'A', 'B', 'C', or 'D' and run the cell." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "You forgot the augmented column. Try again.\n" ] } ], "source": [ "print(quiz3.matrix_quiz(answer = 'B'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Quick Review\n", "\n", "[lab3:sec:quick]: (#Quick-Review)\n", "\n", "Here is a review on the basic matrix operations, including addition,\n", "subtraction and multiplication. These are important in solving linear\n", "systems.\n", "\n", "Here is a short exercise to see how much you remember:\n", "\n", "Let $x = \\left[ \\begin{array}{r} 2 \\\\ 2 \\\\ 7 \\end{array} \\right] , \n", " y = \\left[ \\begin{array}{r} -5 \\\\ 1 \\\\ 3 \\end{array} \\right] , \n", " A = \\left[ \\begin{array}{rrr} 3 & -2 & 10 \\\\ \n", " -6 & 7 & -4 \n", " \\end{array} \\right],\n", " B = \\left[ \\begin{array}{rr} -6 & 4 \\\\ 7 & -1 \\\\ 2 & 9 \n", " \\end{array} \\right]$.\n", "\n", "Calculate the following:\n", "\n", "1. $x + y$\n", "\n", "2. $x^{T}y$\n", "\n", "3. $y-x$\n", "\n", "4. $Ax$\n", "\n", "5. $y^{T}A$\n", "\n", "6. $AB$\n", "\n", "7. $BA$\n", "\n", "8. $AA$\n", "\n", "The solutions to these exercises are available [here](lab3_files/quizzes/quick/quick.html)\n", "\n", "After solving the questions by hand, you can also use Python to check\n", "your answers.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gaussian Elimination\n", "\n", "[lab3:sec:gaus]: <#Gaussian-Elimination> \"Gaussian Elimination\"\n", "\n", "The simplest method for solving a linear system is Gaussian elimination,\n", "which uses three types of *elementary row operations*:\n", "\n", "- Multiplying a row by a non-zero constant ($kE_{ij}$)\n", "\n", "- Adding a multiple of one row to another ($E_{ij} + kE_{kj}$)\n", "\n", "- Exchanging two rows ($E_{ij} \\leftrightarrow E_{kj}$)\n", "\n", "Each row operation corresponds to a step in the solution of the [System\n", "of Equations](#lab3:eq:system) where the equations are combined\n", "together. *It is important to note that none of those operations changes\n", "the solution.* There are two parts to this method: elimination and\n", "back-substitution. The purpose of the process of elimination is to\n", "eliminate the matrix entries below the main diagonal, using row\n", "operations, to obtain a upper triangular matrix with the augmented\n", "column. Then, you will be able to proceed with back-substitution to find\n", "the values of the unknowns.\n", "\n", "Try to solve this set of linear equations:\n", "\n", "$$\\begin{array}{lrcrcrcr}\n", " E_{1j}: & 2x_{1} & + & 8x_{2} & - & 5x_{3} & = & 53 \\\\\n", " E_{2j}: & 3x_{1} & - & 6x_{2} & + & 4x_{3} & = & -48 \\\\\n", " E_{3j}: & x_{1} & + & 2x_{2} & - & x_{3} & = & 13\n", "\\end{array}$$\n", "\n", "The solution to this problem is available [here](lab3_files/quizzes/gaus/gaus.html)\n", "\n", "After solving the system by hand, you can use Python to check your\n", "answer.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Decomposition \n", "[lab3:sec:decomp]: <#Decomposition> \"Decomposition\"\n", "\n", "Any invertible, square matrix, $A$, can be factored out into a product\n", "of a lower and an upper triangular matrices, $L$ and $U$, respectively,\n", "so that $A$ = $LU$. The $LU$- *decomposition* is closely linked to the\n", "process of Gaussian elimination.\n", "\n", "#### Example One\n", "\n", "> Using the matrix from the system of the previous section (Sec\n", "[Gaussian Elimination](#Gaussian-Elimination)), we have:\n", "\n", "> $$A = \\left[ \\begin{array}{rrr} 2 & 8 & -5 \\\\\n", " 3 & -6 & 4 \\\\\n", " 1 & 2 & -1 \\end{array} \\right]$$\n", "\n", "> The upper triangular matrix $U$ can easily be calculated by applying\n", "Gaussian elimination to $A$:\n", "\n", "> $$\\begin{array}{cl}\n", " \\begin{array}{c} E_{2j}-\\frac{3}{2}E_{1j} \\\\\n", " E_{3j}-\\frac{1}{2}E_{1j} \\\\\n", " \\rightarrow \\end{array}\n", " & \n", " \\left[ \\begin{array}{rrr} 2 & 8 & -5 \\\\\n", " 0 & -18 & \\frac{23}{2} \\\\\n", " 0 & -2 & \\frac{3}{2} \n", " \\end{array} \\right] \\\\ \\\\\n", " \\begin{array}{c} E_{3j}-\\frac{1}{9}E_{2j} \\\\\n", " \\rightarrow \\end{array}\n", " &\n", " \\left[ \\begin{array}{rrr} 2 & 8 & -5 \\\\\n", " 0 & -18 & \\frac{23}{2} \\\\\n", " 0 & 0 & \\frac{2}{9}\n", " \\end{array} \\right] = U\n", "\\end{array}$$\n", "\n", "> Note that there is no row exchange.\n", "\n", "> The lower triangular matrix $L$ is calculated with the steps which lead\n", "us from the original matrix to the upper triangular matrix, i.e.:\n", "\n", "> $$\\begin{array}{c} E_{2j}-\\frac{3}{2}E_{1j} \\\\ \\\\\n", " E_{3j}-\\frac{1}{2}E_{1j} \\\\ \\\\\n", " E_{3j}-\\frac{1}{9}E_{2j} \n", "\\end{array}$$\n", "\n", "> Note that each step is a multiple $\\ell$ of equation $m$ subtracted from\n", "equation $n$. Each of these steps, in fact, can be represented by an\n", "elementary matrix. $U$ can be obtained by multiplying $A$ by this\n", "sequence of elementary matrices.\n", "\n", "> Each of the elementary matrices is composed of an identity matrix the\n", "size of $A$ with $-\\ell$ in the ($m,n$) entry. So the steps become:\n", "\n", "> $$\\begin{array}{ccc} \n", " E_{2j}-\\frac{3}{2}E_{1j} & \n", " \\rightarrow &\n", " \\left[ \\begin{array}{rcc} 1 & 0 & 0 \\\\\n", " -\\frac{3}{2} & 1 & 0 \\\\\n", " 0 & 0 & 1\n", " \\end{array} \\right] = R \\\\ \\\\\n", " E_{3j}-\\frac{1}{2}E_{1j} &\n", " \\rightarrow &\n", " \\left[ \\begin{array}{rcc} 1 & 0 & 0 \\\\\n", " 0 & 1 & 0 \\\\\n", " -\\frac{1}{2} & 0 & 1 \n", " \\end{array} \\right] = S \\\\ \\\\\n", " E_{3j}-\\frac{1}{9}E_{2j} &\n", " \\rightarrow &\n", " \\left[ \\begin{array}{crc} 1 & 0 & 0 \\\\\n", " 0 & 1 & 0 \\\\\n", " 0 & -\\frac{1}{9} & 1\n", " \\end{array} \\right] = T\n", "\\end{array}$$\n", "\n", "> and $TSRA$ = $U$. Check this with Python.\n", "\n", "> To get back from $U$ to $A$, the inverse of $R$, $S$ and $T$ are\n", "multiplied onto $U$:\n", "\n", "> $$\\begin{array}{rcl}\n", " T^{-1}TSRA & = & T^{-1}U \\\\\n", " S^{-1}SRA & = & S^{-1}T^{-1}U \\\\ \n", " R^{-1}RA & = & R^{-1}S^{-1}T^{-1}U \\\\\n", "\\end{array}$$\n", "\n", "> So $A$ = $R^{-1}S^{-1}T^{-1}U$. Recall that $A$ = $LU$. If\n", "$R^{-1}S^{-1}T^{-1}$ is a lower triangular matrix, then it is $L$.\n", "\n", "> The inverse of the elementary matrix is the same matrix with only one\n", "difference, and that is, $\\ell$ is in the $a_{mn}$ entry instead of\n", "$-\\ell$. So:\n", "\n", "> $$\\begin{array}{rcl}\n", " R^{-1} & = & \\left[ \\begin{array}{rrr} 1 & 0 & 0 \\\\\n", " \\frac{3}{2} & 1 & 0 \\\\\n", " 0 & 0 & 1\n", " \\end{array} \\right] \\\\ \\\\\n", " S^{-1} & = & \\left[ \\begin{array}{rrr} 1 & 0 & 0 \\\\\n", " 0 & 1 & 0 \\\\\n", " \\frac{1}{2} & 0 & 1\n", " \\end{array} \\right] \\\\ \\\\\n", " T^{-1} & = & \\left[ \\begin{array}{rrr} 1 & 0 & 0 \\\\\n", " 0 & 1 & 0 \\\\\n", " 0 & \\frac{1}{9} & 1\n", " \\end{array} \\right] \n", "\\end{array}$$\n", "\n", "> Multiplying $R^{-1}S^{-1}T^{-1}$ together, we have:\n", "\n", "> $$\\begin{array}{rcl} R^{-1}S^{-1}T^{-1} \n", " & = &\n", " \\left[ \\begin{array}{ccc} 1 & 0 & 0 \\\\\n", " \\frac{3}{2} & 1 & 0 \\\\\n", " \\frac{1}{2} & \\frac{1}{9} & 1\n", " \\end{array} \\right] = L\n", "\\end{array}$$\n", "\n", "> So $A$ is factored into two matrices $L$ and $U$, where\n", "\n", "> $$\\begin{array}{ccc}\n", " L = \\left[ \\begin{array}{ccc} 1 & 0 & 0 \\\\\n", " \\frac{3}{2} & 1 & 0 \\\\\n", " \\frac{1}{2} & \\frac{1}{9} & 1\n", " \\end{array} \\right]\n", "& \\mbox{ and } &\n", " U = \\left[ \\begin{array}{ccc} 2 & 8 & -5 \\\\\n", " 0 & -18 & \\frac{23}{2} \\\\\n", " 0 & 0 & \\frac{2}{9}\n", " \\end{array} \\right]\n", "\\end{array}$$\n", "\n", "> Use Python to confirm that $LU$ = $A$.\n", "\n", "The reason decomposition is introduced here is not because of Gaussian\n", "elimination $-$ one seldom explicitly computes the $LU$ decomposition of\n", "a matrix. However, the idea of factoring a matrix is important for other\n", "direct methods of solving linear systems (of which Gaussian elimination\n", "is only one) and for methods for finding eigenvalues ([Characteristic Equation](#Characteristic-Equation)).\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Round-off Error\n", "\n", "[lab3:sec:round-off-error]: <#Round-off-Error> \"Round-off-Error\"\n", "\n", "When a number is represented in its floating point form, i.e. an\n", "approximation of the number, the resulting error is the *round-off\n", "error*. The floating-point representation of numbers  and the consequent\n", "effects of round-off error were discussed already in Lab \\#2.\n", "\n", "When round-off errors are present in the matrix $A$ or the right hand\n", "side $b$, the linear system $Ax = b$ may or may not give a solution that\n", "is close to the real answer. When a matrix $A$ “magnifies” the effects\n", "of round-off errors in this way, we say that $A$ is an ill-conditioned\n", "matrix.\n", "\n", "#### Example Two\n", "[lab3:eg:round]: <#Example-Two> \"Example Two\"\n", "\n", "> Let’s see an example:\n", "\n", "> Suppose\n", "\n", "> $$A = \\left[ \\begin{array}{cc} 1 & 1 \\\\ 1 & 1.0001 \n", " \\end{array} \\right]$$\n", "\n", "> and consider the system:\n", "\n", ">
\n", "(Ill-conditioned version one):\n", " $$\\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 1 \\\\ 1 & 1.0001 \\end{array} \n", "&\n", " \\left| \\begin{array}{c} 2 \\\\ 2 \\end{array} \\right]\n", "\\end{array} \\right.\n", "$$\n", "
\n", "\n", "> The condition number, $K$, of a matrix, defined in [Condition Number](#Condition-Number), is a measure of how well-conditioned a matrix is. If\n", "$K$ is large, then the matrix is ill-conditioned, and Gaussian\n", "elimination will magnify the round-off errors. The condition number of\n", "$A$ is 40002. You can use Python to check this number.\n", "\n", "> The solution to this is $x_1$ = 2 and $x_2$ = 0. However, if the system\n", "is altered a little as follows:\n", "\n", ">
\n", "(Ill-conditioned version two):\n", " $$\\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 1 \\\\ 1 & 1.0001 \\end{array}\n", "&\n", " \\left| \\begin{array}{c} 2 \\\\ 2.0001 \\end{array} \\right]\n", "\\end{array} \\right.\n", "$$\n", "
\n", "\n", "> Then, the solution becomes $x_1$ = 1 and $x_2$ = 1. A change in the\n", "fifth significant digit was amplified to the point where the solution is\n", "not even accurate to the first significant digit. $A$ is an\n", "ill-conditioned matrix. You can set up the systems [Ill-conditioned version one](#lab3:eq:illbefore)\n", "and [Ill-conditioned version two](#lab3:eq:illafter) in Python, and check the answers yourself." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example Three\n", "\n", "[lab3:eg:inacc]: <#Example-Three> \"Example Three\"\n", "\n", "> Use Python to try the following example. First solve the\n", "system $A^{\\prime}x = b$; then solve $A^{\\prime}x = b2$. Find the\n", "condition number of $A^{\\prime}$.\n", "\n", ">$$\\begin{array}{ccccc}\n", "A^{\\prime} = \\left[ \\begin{array}{cc} 0.0001 & 1 \\\\ 1 & 1 \n", " \\end{array} \\right]\n", "& , & \n", "b = \\left[ \\begin{array}{c} 1 \\\\ 2 \\end{array} \\right]\n", "& \\mbox{and} &\n", "b2 = \\left[ \\begin{array}{c} 1 \\\\ 2.0001 \\end{array} \\right] .\n", "\\end{array}$$\n", "\n", "> You will find that the solution for $A^{\\prime}x = b$ is $x_1$ = 1.0001\n", "and $x_2$ = 0.9999, and the solution for $A^{\\prime}x = b2$ is $x_1$ =\n", "1.0002 and $x_2$ = 0.9999 . So a change in $b$ did not result in a large\n", "change in the solution. Therefore, $A^{\\prime}$ is a well-conditioned\n", "matrix. In fact, the condition number is approximately 2.6.\n", "\n", "> Nevertheless, even a well conditioned system like $A^{\\prime}x =\n", "b$ leads to inaccuracy if the wrong solution method is used, that is, an\n", "algorithm which is sensitive to round-off error. If you use Gaussian\n", "elimination to solve this system, you might be misled that $A^{\\prime}$\n", "is ill-conditioned. Using Gaussian elimination to solve $A^{\\prime}x=b$:\n", "\n", "> $$\\begin{array}{cl}\n", "& \\left[\n", "\\begin{array}{cc}\n", " \\begin{array}{cc} 0.0001 & 1 \\\\ 1 & 1 \\end{array} \n", "&\n", " \\left|\n", " \\begin{array}{c} 1 \\\\ 2 \\end{array} \\right]\n", "\\end{array} \\right. \\\\ \\\\\n", "\\begin{array}{c} 10,000E_{1j} \\\\ \\rightarrow \\end{array} &\n", "\\left[\n", "\\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 10,000 \\\\ 1 & 1 \\end{array} \n", "&\n", " \\left|\n", " \\begin{array}{c} 10,000 \\\\ 2 \\end{array} \\right]\n", "\\end{array} \\right. \\\\ \\\\\n", "\\begin{array}{c} E_{2j}-E_{1j} \\\\ \\rightarrow \\end{array} &\n", "\\left[\n", "\\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 10,000 \\\\ 0 & -9,999 \\end{array} \n", "&\n", " \\left|\n", " \\begin{array}{c} 10,000 \\\\ -9,998 \\end{array} \\right]\n", "\\end{array} \\right.\n", "\\end{array}$$\n", "\n", "> At this point, if you continue to solve the system as is, you will get\n", "the expected answers. You can check this with Python. However, if you\n", "make changes to the matrix here by rounding -9,999 and -9,998 to -10,000, \n", "the final answers will be different:\n", "\n", "> $$\\begin{array}{cl}\n", "& \\left[\n", "\\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 10,000 \\\\ 0 & -10,000 \\end{array}\n", "&\n", " \\left|\n", " \\begin{array}{c} 10,000 \\\\ -10,000 \\end{array} \\right]\n", "\\end{array} \\right. \n", "\\end{array}$$\n", "\n", "> The result is $x_1$ = 0 and $x_2$ = 1, which is quite different from the\n", "correct answers. So Gaussian elimination might mislead you to think that\n", "a matrix is ill-conditioned by giving an inaccurate solution to the\n", "system. In fact, the problem is that Gaussian elimination on its own is\n", "a method that is unstable in the presence of round-off error, even for\n", "well-conditioned matrices. Can this be fixed?\n", "\n", "> You can try the example with Python.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partial Pivoting\n", "\n", "There are a number of ways to avoid inaccuracy, one of which is applying\n", "partial pivoting to the Gaussian elimination.\n", "\n", "\n", "Consider the example from the previous section. In order to avoid\n", "multiplying by 10,000, another pivot is desired in place of 0.0001. The\n", "goal is to examine all the entries in the first column, find the entry\n", "that has the largest value, and exchange the first row with the row that\n", "contains this element. So this entry becomes the pivot. This is partial\n", "pivoting. Keep in mind that switching rows is an elementary operation\n", "and has no effect on the solution.\n", "\n", "\n", "In the original Gaussian elimination algorithm, row exchange is done\n", "only if the pivot is zero. In partial pivoting, row exchange is done so\n", "that the largest value in a certain column is the pivot. This helps to\n", "reduce the amplification of round-off error." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example Four\n", "\n", "> In the matrix $A^{\\prime}$ from [Example Two](#Example-Two), 0.0001 from\n", "column one is the first pivot. Looking at this column, the entry, 1, in\n", "the second row is the only other choice in this column. Obviously, 1 is\n", "greater than 0.0001. So the two rows are exchanged.\n", "\n", "> $$\\begin{array}{cl}\n", " & \\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 0.0001 & 1 \\\\ 1 & 1 \\end{array} \n", " & \\left|\n", " \\begin{array}{c} 1 \\\\ 2 \\end{array} \\right]\n", " \\end{array} \\right. \\\\ \\\\\n", " \\begin{array}{c} E_{1j} \\leftrightarrow E_{2j} \\\\ \\rightarrow\n", " \\end{array}\n", " & \\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 1 \\\\ 0.0001 & 1 \\end{array} \n", " & \\left|\n", " \\begin{array}{c} 2 \\\\ 1 \\end{array} \\right] \n", " \\end{array} \\right. \\\\ \\\\\n", " \\begin{array}{c} E_{2j}-0.0001E_{1j} \\\\ \\rightarrow \\end{array} \n", " & \\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 1 \\\\ 0 & 0.9999 \\end{array} \n", " & \\left|\n", " \\begin{array}{c} 2 \\\\ 0.9998 \\end{array} \\right]\n", " \\end{array} \\right.\n", "\\end{array}$$\n", "\n", "> The same entries are rounded off:\n", "\n", "> $$\\begin{array}{cl}\n", " & \\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 1 \\\\ 0 & 1 \\end{array}\n", " &\n", " \\left|\n", " \\begin{array}{c} 2 \\\\ 1 \\end{array} \\right]\n", " \\end{array} \\right. \\\\ \\\\\n", " \\begin{array}{c} E_{1j}-E_{2j} \\\\ \\rightarrow \\end{array}\n", " & \\left[ \\begin{array}{cc}\n", " \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array}\n", " &\n", " \\left|\n", " \\begin{array}{c} 1 \\\\ 1 \\end{array} \\right]\n", " \\end{array} \\right.\n", "\\end{array}$$\n", "\n", "> So the solution is $x_1$ = 1 and $x_2$ = 1, and this is a close\n", "approximation to the original solution, $x_1$ = 1.0001 and $x_2$ =\n", "0.9999.\n", "\n", "\n", "> You can try the example with Python.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n", "\n", "\n", "\n", "Note: This section has described row pivoting. The same process can be\n", "applied to columns, with the resulting procedure being called column\n", "pivoting." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Full Pivoting\n", "\n", "Another way to get around inaccuracy (Example [Three](#Example-Three)) is to\n", "use Gaussian elimination with full pivoting. Sometimes, even partial\n", "pivoting can lead to problems. With full pivoting, in addition to row\n", "exchange, columns will be exchanged as well. The purpose is to use the\n", "largest entries in the whole matrix as the pivots.\n", "\n", "#### Example Five\n", "\n", "> Given the following:\n", "\n", "> $$ \\begin{array}{cccc}\n", "& A^{''} = \\left[ \\begin{array}{ccc} 0.0001 & 0.0001 & 0.5 \\\\\n", " 0.5 & 1 & 1 \\\\\n", " 0.0001 & 1 & 0.0001 \n", " \\end{array} \n", " \\right] \n", "& \\ \\ \\ &\n", "b^{'} = \\left[ \\begin{array}{c} 1 \\\\ 0 \\\\ 1 \n", " \\end{array} \n", " \\right] \\\\ \\\\\n", "\\rightarrow &\n", "\\left[ \\begin{array}{cc} \n", " \\begin{array}{ccc} 0.0001 & 0.0001 & 0.5 \\\\\n", " 0.5 & 1 & 1 \\\\\n", " 0.0001 & 1 & 0.0001\n", " \\end{array}\n", " & \\left| \\begin{array}{c} 1 \\\\ 0 \\\\ 1 \n", " \\end{array} \n", " \\right]\n", " \\end{array} \n", "\\right. & & \n", "\\end{array}$$\n", "\n", "> Use Python to find the condition number of $A^{''}$ and the solution to\n", "this system.\n", "\n", "> Looking at the system, if no rows are exchanged, then taking 0.0001 as\n", "the pivot will magnify any errors made in the elements in column 1 by a\n", "factor of 10,000. With partial pivoting, the first two rows can be\n", "exchanged (as below):\n", "\n", "> $$\\begin{array}{cl}\n", "\\begin{array}{c} \\\\ E_{1j} \\leftrightarrow E_{2j} \\\\ \\rightarrow\n", "\\end{array} \n", "& \\begin{array}{ccccc}\n", " & x_1 & x_2 & x_3 & \\\\\n", " \\left[ \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right. \n", " & \\begin{array}{c} 0.5 \\\\ 0.0001 \\\\ 0.0001 \\end{array} \n", " & \\begin{array}{c} 1 \\\\ 0.0001 \\\\ 1 \\end{array}\n", " & \\begin{array}{c} 1 \\\\ 0.5 \\\\ 0.0001 \\end{array} \n", " & \\left| \\begin{array}{c} 0 \\\\ 1 \\\\ 1 \\end{array} \\right]\n", " \\end{array}\n", "\\end{array} $$ \n", "\n", "\n", "> and the magnification by 10,000 is avoided. Now the matrix will be\n", "expanded by a factor of 2. However, if the entry 1 is used as the pivot,\n", "then the matrix does not need to be expanded by 2 either. The only way\n", "to put 1 in the position of the first pivot is to perform a column\n", "exchange between columns one and two, or between columns one and three.\n", "This is full pivoting.\n", "\n", "> Note that when columns are exchanged, the variables represented by the\n", "columns are switched as well, i.e. when columns one and two are\n", "exchanged, the new column one represents $x_2$ and the new column two\n", "represents $x_1$. So, we must keep track of the columns when performing\n", "column pivoting.\n", "\n", ">So the columns one and two are exchanged, and the matrix becomes:\n", "\n", "> $$\\begin{array}{cl}\n", "\\begin{array}{c} \\\\ E_{i1} \\leftrightarrow E_{i2} \\\\ \\rightarrow\n", "\\end{array}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_1 & x_3 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0.0001 \\\\ 1 \\end{array}\n", " & \\begin{array}{c} 0.5 \\\\ 0.0001 \\\\ 0.0001 \\end{array} \n", " & \\begin{array}{c} 1 \\\\ 0.5 \\\\ 0.0001 \\end{array} \n", " & \\left| \\begin{array}{c} 0 \\\\ 1 \\\\ 1 \\end{array} \\right]\n", "\\end{array} \\\\ \\\\\n", "\\begin{array}{c} \\\\ E_{2j}-0.0001E_{1j} \\\\ E_{3j}-E_{1j} \\\\ \\rightarrow\n", "\\end{array}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_1 & x_3 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0.5 \\\\ 0.00005 \\\\ -0.4999 \\end{array}\n", " & \\begin{array}{c} 1 \\\\ 0.4999 \\\\ -0.9999 \\end{array}\n", " & \\left| \\begin{array}{c} 0 \\\\ 1 \\\\ 1 \\end{array} \\right]\n", "\\end{array}\n", "\\end{array}$$\n", "\n", "> If we assume rounding is performed, then the entries are rounded off:\n", "\n", "> $$\\begin{array}{cl}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_1 & x_3 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0.5 \\\\ 0.00005 \\\\ -0.5 \\end{array}\n", " & \\begin{array}{c} 1 \\\\ 0.5 \\\\ -1 \\end{array}\n", " & \\left| \\begin{array}{c} 0 \\\\ 1 \\\\ 1 \\end{array} \\right]\n", "\\end{array} \\\\ \\\\\n", "\\begin{array}{c} \\\\ -E_{3j} \\\\\n", " E_{2j} \\leftrightarrow E_{3j} \\\\ \n", " E_{i2} \\leftrightarrow E_{i3} \\\\ \n", " \\rightarrow\n", "\\end{array}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_3 & x_1 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 1 \\\\ 1 \\\\ 0.5 \\end{array}\n", " & \\begin{array}{c} 0.5 \\\\ 0.5 \\\\ 0.00005 \\end{array}\n", " & \\left| \\begin{array}{r} 0 \\\\ -1 \\\\ 1 \\end{array} \\right]\n", "\\end{array} \\\\ \\\\\n", "\\begin{array}{c} E_{1j}-E_{2j} \\\\ \n", " E_{3j}-0.5 E_{2j} \\\\\n", " \\rightarrow\n", "\\end{array}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_3 & x_1 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0 \\\\ 1 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0 \\\\ 0.5 \\\\ -0.24995 \\end{array}\n", " & \\left| \\begin{array}{r} 1 \\\\ -1 \\\\ 1.5 \\end{array} \\right]\n", "\\end{array}\n", "\\end{array}$$\n", "\n", "> Rounding off the matrix again:\n", "\n", "> $$\\begin{array}{cl}\n", "\\begin{array}{c} \\rightarrow \\end{array}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_3 & x_1 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0 \\\\ 1 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0 \\\\ 0.5 \\\\ -0.25 \\end{array}\n", " & \\left| \\begin{array}{r} 1 \\\\ -1 \\\\ 1.5 \\end{array} \\right]\n", "\\end{array} \\\\ \\\\ \n", "\\begin{array}{c} E_{2j}-2E_{3j} \\\\ \n", " 4E_{3j} \\\\\n", " \\rightarrow \\end{array}\n", "& \\begin{array}{ccccc}\n", " & x_2 & x_3 & x_1 & \\\\\n", " \\left. \\begin{array}{c} \\\\ \\\\ \\\\ \\end{array} \\right[\n", " & \\begin{array}{c} 1 \\\\ 0 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0 \\\\ 1 \\\\ 0 \\end{array}\n", " & \\begin{array}{c} 0 \\\\ 0 \\\\ 1 \\end{array}\n", " & \\left| \\begin{array}{r} 1 \\\\ 2 \\\\ -6 \\end{array} \\right]\n", "\\end{array}\n", "\\end{array}$$\n", "\n", "> So reading from the matrix, $x_1$ = -6, $x_2$ = 1 and $x_3$ = 2. Compare\n", "this with the answer you get with Python, which is $x_1\n", "\\approx$ -6.0028, $x_2 \\approx$ 1.0004 and $x_3 \\approx$ 2.0010 .\n", "\n", "> Using full pivoting with Gaussian elimination, expansion of the error by\n", "large factors is avoided. In addition, the approximated solution, using\n", "rounding (which is analogous to the use of floating point\n", "approximations), is close to the correct answer.\n", "\n", "> You can try the example with Python.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Summary\n", "\n", "In a system $Ax = b$, if round-off errors in $A$ or $b$ affect the\n", "system such that it may not give a solution that is close to the real\n", "answer, then $A$ is ill-conditioned, and it has a very large condition\n", "number.\n", "\n", "Sometimes, due to a poor algorithm, such as Gaussian elimination without\n", "pivoting, a matrix may appear to be ill-conditioned, even though it is\n", "not. By applying partial pivoting, this problem is reduced, but partial\n", "pivoting will not always eliminate the effects of round-off error. An\n", "even better way is to apply full pivoting. Of course, the drawback of\n", "this method is that the computation is more expensive than plain\n", "Gaussian elimination.\n", "\n", "An important point to remember is that partial and full pivoting\n", "minimize the effects of round-off error for well-conditioned matrices.\n", "If a matrix is ill-conditioned, these methods will not provide a\n", "solution that approximates the real answer. As an exercise, you can try\n", "to apply full pivoting to the ill-conditioned matrix $A$ seen at the\n", "beginning of this section (Example [Two](#Example-Two)). You will find that\n", "the solution is still inaccurate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrix Inversion\n", "\n", "Given a square matrix $A$. If there is a matrix that will cancel $A$,\n", "then it is the *inverse* of $A$. In other words, the matrix, multiplied\n", "by its inverse, will give the identity matrix $I$.\n", "\n", "Try to find the inverse for the following matrices:\n", "\n", "1. $$A = \\left[ \\begin{array}{rrr} 1 & -2 & 1 \\\\\n", " 3 & 1 & -1 \\\\\n", " -1 & 9 & -5 \\end{array} \\right]$$\n", "\n", "2. $$B = \\left[ \\begin{array}{rrr} 5 & -2 & 4 \\\\\n", " -3 & 1 & -5 \\\\\n", " 2 & -1 & 3 \\end{array} \\right]$$\n", "\n", "The solutions to these exercises are available [here](lab3_files/quizzes/inverse/inverse.html)\n", "\n", "After solving the questions by hand, you can use Python to check your\n", "answer.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Determinant\n", "\n", "Every square matrix $A$ has a scalar associated with it. This number is\n", "the determinant of the matrix, represented as $\\det(A)$. Its absolute\n", "value is the volume of the parallelogram that can be generated from the\n", "rows of $A$.\n", "\n", "A few special properties to remember about determinants:\n", "\n", "1. $A$ must be a square matrix.\n", "\n", "2. If $A$ is singular, $\\det(A) =0$, i.e. $A$ does not have an inverse.\n", "\n", "3. The determinant of a $2 \\times 2$ matrix is just the difference\n", " between the products of the diagonals, i.e.\n", "\n", " $$\\begin{array}{ccc}\n", " \\left[ \\begin{array}{cc} a & b \\\\ c & d \\end{array} \\right] &\n", " = &\n", " \\begin{array}{ccc} ad & - & bc \\end{array}\n", " \\end{array}$$\n", "\n", "4. For any diagonal, upper triangular or lower triangular matrix $A$,\n", " $\\det(A)$ is the product of all the entries on the diagonal,\n", " \n", "#### Example Six\n", "\n", "> $$\\begin{array}{cl}\n", " & \\det \\left[\n", " \\begin{array}{rrr}\n", " 2 & 5 & -8 \\\\\n", " 0 & 1 & 7 \\\\\n", " 0 & 0 & -4 \n", " \\end{array} \\right] \\\\ \\\\\n", " = & 2 \\times 1 \\times -4 \\\\\n", " = & -8\n", " \\end{array}$$ \n", " \n", "> Graphically, the parallelogram looks as follows:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The basic procedure in finding a determinant of a matrix larger than 2\n", "$\\times$ 2 is to calculate the product of the non-zero entries in each\n", "of the *permutations* of the matrix, and then add them together. A\n", "permutation of a matrix $A$ is a matrix of the same size with one\n", "element from each row and column of $A$. The sign of each permutation is\n", "$+$ or $-$ depending on whether the permutation is odd or even. This is\n", "illustrated in the following example …" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example Seven\n", "\n", "> $$I = \\left[ \\begin{array}{ccc} 1 & 2 & 3 \\\\\n", " 4 & 5 & 6 \\\\\n", " 7 & 8 & 9 \\end{array} \\right]$$\n", "\n", "> will have the following permutations:\n", "\n", "> $$\\begin{array}{cccccc}\n", "+\\left[ \\begin{array}{ccc} 1 & & \\\\\n", " & 5 & \\\\\n", " & & 9 \\end{array} \\right]\n", "& , &\n", "+\\left[ \\begin{array}{ccc} & 2 & \\\\\n", " & & 6 \\\\\n", " 7 & & \\end{array} \\right]\n", "& , &\n", "+\\left[ \\begin{array}{ccc} & & 3 \\\\\n", " 4 & & \\\\\n", " & 8 & \\end{array} \\right]\n", "& , \\\\ \\\\\n", "-\\left[ \\begin{array}{ccc} 1 & & \\\\\n", " & & 6 \\\\\n", " & 8 & \\end{array} \\right]\n", "& , &\n", "-\\left[ \\begin{array}{ccc} & 2 & \\\\\n", " 4 & & \\\\\n", " & & 9 \\end{array} \\right]\n", "& , &\n", "-\\left[ \\begin{array}{ccc} & & 3 \\\\\n", " & 5 & \\\\\n", " 7 & & \\end{array} \\right]\n", "& .\n", "\\end{array}$$\n", "\n", "> The determinant of the above matrix is then given by \n", "$$\\begin{array}{ccl}\n", " det(A) & = & +1\\cdot 5\\cdot 9 + 2 \\cdot 6 \\cdot 7 + 3 \\cdot 4 \\cdot 8 - 1\n", "\\cdot 6 \\cdot 8 - 2 \\cdot 4 \\cdot 9 - 3 \\cdot 5 \\cdot 7 \\\\\n", "& = & 0\\end{array}$$\n", "\n", "For each of the following matrices, determine whether or not it has an\n", "inverse:\n", "\n", "1. $$A = \\left[ \\begin{array}{rrr} 3 & -2 & 1 \\\\\n", " 1 & 5 & -1 \\\\\n", " -1 & 0 & 0 \\end{array} \\right]$$\n", "\n", "2. $$B = \\left[ \\begin{array}{rrr} 4 & -6 & 1 \\\\\n", " 1 & -3 & 1 \\\\\n", " 2 & 0 & -1 \\end{array} \\right]$$\n", "\n", "3. Try to solve this by yourself first, and use Python to check your\n", " answer:\n", "\n", " $$C = \\left[ \\begin{array}{rrrr} 4 & -2 & -7 & 6 \\\\\n", " -3 & 0 & 1 & 0 \\\\\n", " -1 & -1 & 5 & -1 \\\\\n", " 0 & 1 & -5 & 3\n", " \\end{array} \\right]$$\n", "\n", "The solutions to these exercises are available [here](lab3_files/quizzes/det/det.html)\n", "\n", "After solving the questions by hand, you can use Python to check your\n", "answer.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computational cost of Gaussian elimination\n", "\n", "[lab3:sec:cost]: <#Computational-cost-of-Gaussian-elimination> \"Computational Cost\"\n", "\n", "Although Gaussian elimination is a basic and relatively simple technique\n", "to find the solution of a linear system, it is a costly algorithm. Here\n", "is an operation count of this method:\n", "\n", "For a $n \\times n$ matrix, there are two kinds of operations to\n", "consider:\n", "\n", "1. division ( *div*) - to find the multiplier from a\n", " chosen pivot\n", "\n", "2. multiplication-subtraction ( *mult/sub* ) - to\n", " calculate new entries for the matrix\n", "\n", "Note that an addition or subtraction operation has negligible cost in\n", "relation to a multiplication or division operation. So the subtraction\n", "in this case can be treated as one with the multiplication operation.\n", "The first pivot is selected from the first row in the matrix. For each\n", "of the remaining rows, one div and $(n-1)$ mult/sub operations are used\n", "to find the new entries. So there are $n$ operations performed on each\n", "row. With $(n-1)$ rows, there are a total of $n(n-1) = n^{2}-n$\n", "operations associated with this pivot.\n", "\n", "Since the subtraction operation has negligible cost in relation to the\n", "multiplication operation, there are $(n-1)$ operations instead of\n", "$2(n-1)$ operations.\n", "For the second pivot, which is selected from the second row of the\n", "matrix, similar analysis is applied. With the remaining $(n-1)\n", "\\times (n-1)$ matrix, each row has one div and $(n-2)$ mult/sub\n", "operations. For the whole process, there are a total of $(n-1)(n-2) =\n", "(n-1)^{2} - (n-1)$ operations.\\\n", "\n", "For the rest of the pivots, the number of operations for a remaining\n", "$k \\times k$ matrix is $k^{2} - k$.\\\n", "\n", "The following is obtained when all the operations are added up:\n", "\n", "$$\\begin{array}{l} \n", "(1^{2}+\\ldots +n^{2}) - (1+\\ldots +n) \\\\ \\\\\n", "= \\frac{n(n+1)(2n+1)}{6} - \\frac{n(n+1)}{2} \\\\ \\\\\n", "= \\frac{n^{3}-n}{3} \\\\ \\\\\n", "\\approx O(n^{3}) \n", "\\end{array}$$\n", "\n", "As one can see, the Gaussian elimination is an $O(n^{3})$ algorithm. For\n", "large matrices, this can be prohibitively expensive. There are other\n", "methods which are more efficient, e.g. see [Iterative Methods](#Iterative-Methods). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem One\n", "\n", "[lab3:sec2:carbon]:<#Problem-One>\n", "\n", "Consider a very simple three box model of the movement of a pollutant in\n", "the atmosphere, fresh-water and ocean. The mass of the atmosphere is MA\n", "(5600 x 10$^{12}$ tonnes), the mass of the fresh-water is MF (360 x\n", "10$^{12}$tonnes) and the mass of the upper layers of the ocean is MO\n", "(50,000 x 10$^{12}$ tonnes). The amount of pollutant in the atmosphere\n", "is A, the amount in the fresh water is F and the amount in the ocean is\n", "O.\n", "\n", "The pollutant is going directly into the atmosphere at a rate P1 = 1000\n", "tonnes/year and into th efresh-water system at a rate P2 = 2000\n", "tonnes/year. The pollutant diffuses between the atmosphere and ocean at\n", "a rate depending linearly on the difference in concentration with a\n", "diffusion constant L1 = 200 tonnes/year. The diffusion between the\n", "fresh-water system and the atmosphere is faster as the fresh water is\n", "shallower, L2 = 500 tonnes/year. The fresh-water system empties into the\n", "ocean at the rate of Q = 36 x 10$^{12}$ tonnes/year. Lastly the\n", "pollutant decays (like radioactivity) at a rate L3 = 0.05 /year.\n", "\n", "See the graphical presentation of the cycle described above in\n", "Figure [Box Model](#Figure-Box-Model) Schematic for Problem 1.\n", "\n", "- a\\) Consider the steady state. There is no change in A, O, or F. Write\n", " down the three linear governing equations. Write the equations as an\n", " augmented matrix. Use Octave to find the solution.\n", "\n", "- b\\) Show mathematically that there is no solution to this problem with L3\n", " = 0. Why, physically, is there no solution?\n", "\n", "- c\\) Show mathematically that there is an infinite number of solutions if\n", " L3 = 0 and P1 = P2 = 0. Why, physically?\n", "\n", "- d\\) For part c) above, what needs to be specified in order to determine a\n", " single physical solution. How would you put this in the matrix equation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "**Figure Box Model:** Schematic for [Problem One](#Problem-One).\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Eigenvalue Problems\n", "\n", "This section is a review of eigenvalues and eigenvectors.\n", "\n", "### Characteristic Equation\n", "[lab3:sec:eigval]: <#Characteristic-Equation>\n", "\n", "The basic equation for eigenvalue problems is the characteristic\n", "equation, which is:\n", "\n", "$$\\det( A - \\lambda I ) = 0$$\n", "\n", "where $A$ is a square matrix, I is an identity the same size as $A$, and\n", "$\\lambda$ is an eigenvalue of the matrix $A$.\n", "\n", "In order for a number to be an eigenvalue of a matrix, it must satisfy\n", "the characteristic equation, for example:\n", "\n", "#### Example Eight\n", "\n", "> Given\n", "\n", "> $$A = \\left[ \\begin{array}{rr} 3 & -2 \\\\ -4 & 5 \\end{array} \\right]$$\n", "\n", "> To find the eigenvalues of $A$, you need to solve the characteristic\n", "equation for all possible $\\lambda$.\n", "\n", "> $$\\begin{array}{ccl}\n", "0 & = & \\det (A - \\lambda I) \\\\\n", "& = & \\begin{array}{cccc}\n", " \\det & \\left( \n", " \\left[ \\begin{array}{rr} 3 & -2 \\\\ -4 & 5 \\end{array} \\right] \\right. &\n", " - &\n", " \\lambda \\left. \\left[ \\begin{array}{rr} 1 & 0 \\\\ 0 & 1 \\end{array} \\right] \\right)\n", " \\end{array} \\\\ \\\\\n", "& = & \\begin{array}{cc}\n", " \\det & \n", " \\left[ \\begin{array}{cc} 3-\\lambda & -2 \\\\ -4 & 5-\\lambda\n", " \\end{array} \\right]\n", " \\end{array} \\\\ \\\\\n", "& = & \\begin{array}{ccc} (3-\\lambda)(5-\\lambda) & - & (-2)(-4) \n", " \\end{array} \\\\ \\\\\n", "& = & (\\lambda - 1)(\\lambda - 7) \\\\ \\\\\n", "\\end{array}$$\n", "\n", "> So, $\\lambda = 1 \\mbox{ or } 7$, i.e. the eigenvalues of the matrix $A$\n", "are 1 and 7.\n", "\n", "> You can use Python to check this answer.\n", "\n", "Find the eigenvalues of the following matrix:\n", "\n", "$$B = \\left[\n", " \\begin{array}{ccc} 3 & 2 & 4 \\\\ 2 & 0 & 2 \\\\ 4 & 2 & 3\n", " \\end{array} \\right]$$\n", "\n", "The solution to this problem is available [here](lab3_files/quizzes/char/char.html)\n", "\n", "After solving the questions by hand, you can use Python to check your\n", "answer.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Condition Number\n", "[lab3.sec.cond]: \n", "\n", "The eigenvalues of a matrix $A$ can be used to calculate an\n", "approximation to the condition number $K$ of the matrix, i.e.\n", "\n", "$$K \\approx \\left| \\frac{\\lambda_{max}}{\\lambda_{min}} \\right|$$\n", "\n", "where $\\lambda_{max}$ and $\\lambda_{min}$ are the maximum and minimum\n", "eigenvalues of $A$. When $K$ is large, i.e. the $\\lambda_{max}$ and\n", "$\\lambda_{min}$ are far apart, then $A$ is ill-conditioned.\n", "\n", "The mathematical definition of $K$ is\n", "\n", "$$K = \\|A\\|\\|A^{-1}\\|$$\n", "\n", "where $\\|\\cdot\\|$ represents the norm of a matrix.\n", "\n", "There are a few norms which can be chosen for the formula. The default one used\n", "in Python for finding $K$ is the 2-norm of the matrix. To see how to\n", "compute the norm of a matrix, see a linear algebra text. Nevertheless,\n", "the main concern here is the formula, and the fact that this can be very\n", "expensive to compute. Actually, the computing of $A^{-1}$ is the costly\n", "operation.\n", "\n", "Note: In Python, the results from the function *cond*($A$) can have\n", "round-off errors.\n", "\n", "For the matrices in this section (A from Example 8 and B just below it) for which you have \n", "found the\n", "eigenvalues, use the built-in Python function *np.linalg.cond*($A$) to find $K$,\n", "and compare this result with the $K$ approximated from the eigenvalues.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Eigenvectors\n", "\n", "Another way to look at the characteristic equation is using vectors\n", "instead of determinant. For a number to be an eigenvalue of a matrix, it\n", "must satisfy this equation:\n", "\n", "$$( A - \\lambda I ) x = 0$$\n", "\n", "where $A$ is a $n \\times n$ square matrix, $I$ is an identity matrix the\n", "same size as $A$, $\\lambda$ is an eigenvalue of $A$, and $x$ is a\n", "non-zero vector associated with the particular eigenvalue that will make\n", "this equation true. This vector is the eigenvector. The eigenvector is\n", "not necessarily unique for an eigenvalue. This will be further discussed\n", "below after the example.\n", "\n", "The above equation can be rewritten as:\n", "\n", "$$A x = \\lambda x$$\n", "\n", "For each eigenvalue of $A$, there is a corresponding eigenvector. Below\n", "is an example.\n", "\n", "#### Example Nine\n", "\n", "> Following the example from the previous section:\n", "\n", "> $$A = \\left[ \\begin{array}{rr} 3 & -2 \\\\ -4 & 5 \\end{array} \\right]$$\n", "\n", "> The eigenvalues, $\\lambda$, for this matrix are 1 and 7. To find the\n", "eigenvectors for the eigenvalues, you need to solve the equation:\n", "\n", "> $$( A - \\lambda I ) x = 0.$$ \n", "\n", "> This is just a linear system $A^{\\prime}x = b$, where\n", "$A^{\\prime} = ( A - \\lambda I )$, $b = 0$. To find the eigenvectors, you\n", "need to find the solution to this augmented matrix for each $\\lambda$\n", "respectively,\n", "\n", "> $$\\begin{array}{cl}\n", "& ( A - \\lambda I ) x = 0 \\ \\ \\ \\ \\ \\ \\ \\ {\\rm where} \\ \\ \\ \\ \\ \\ \\lambda = 1 \\\\\n", "\\; & \\; \\\\\n", "\\rightarrow & \n", "\\left( \\begin{array}{ccc} \n", " \\left[ \\begin{array}{rr} 3 & -2 \\\\ -4 & 5 \\end{array} \\right] \n", " & - &\n", " 1 \\left[ \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right]\n", "\\end{array} \\right) x = 0 \\\\ \\\\\n", "\\rightarrow &\n", "\\left[ \\begin{array}{cc}\n", " \\begin{array}{rr} 2 & -2 \\\\ -4 & 4 \\end{array}\n", " & \\left| \\begin{array}{c} 0 \\\\ 0 \\end{array} \\right]\n", "\\end{array} \\right. \\\\ \\\\ \n", "\\rightarrow &\n", "\\left[ \\begin{array}{cc}\n", " \\begin{array}{rr} 1 & -1 \\\\ 0 & 0 \\end{array}\n", " & \\left| \\begin{array}{c} 0 \\\\ 0 \\end{array} \\right]\n", "\\end{array} \\right.\n", "\\end{array}$$\n", "\n", "> Reading from the matrix,\n", "\n", "> $$\\begin{array}{ccccc} x_1 & - & x_2 & = & 0 \\\\\n", " & & 0 & = & 0 \\end{array}$$\n", "\n", "> As mentioned before, the eigenvector is not unique for a given\n", "eigenvalue. As seen here, the solution to the matrix is a description of\n", "the direction of the vectors that will satisfy $Ax = \\lambda x$. Letting\n", "$x_1 = 1$, then $x_2 = 1$. So the vector (1, 1) is an eigenvector for\n", "the matrix $A$ when $\\lambda = 1$. (So is (-1,-1), (2, 2), etc)\\\n", "\n", "> In the same way for $\\lambda = 7$, the solution is\n", "\n", "> $$\\begin{array}{ccccc} 2 x_1 & + & x_2 & = & 0 \\\\\n", " & & 0 & = & 0 \\end{array}$$\n", "\n", "> So an eigenvector here is x = (1, -2).\n", "\n", "> Using Python:\n", "\n", " A = np.array([[3, -2], [-4, 5]])\n", " lamb, x = np.linalg.eig(A)\n", " print(lamb)\n", " print (x)\n", "\n", " [ 1. 7.]\n", " [[-0.70710678 0.4472136 ]\n", " [-0.70710678 -0.89442719]]\n", "\n", "> Matrix $x$ is the same size a $A$ and vector $lamb$ is the size of one dimension of $A$. Each column of\n", "$x$ is a unit eigenvector of $A$, and $lamb$ values\n", "are the eigenvalues of $A$. Reading from the result, for $\\lambda$ = 1,\n", "the corresponding unit eigenvector is (-0.70711, -0.70711). The answer\n", "from working out the example by hand is (1, 1), which is a multiple of\n", "the unit eigenvector from Python.\n", "\n", "> (The unit eigenvector is found by dividing the eigenvector by its\n", "magnitude. In this case, $\\mid$(1,1)$\\mid$ = $\\sqrt{1^2 +\n", " 1^2}$ = $\\sqrt{2}$, and so the unit eigenvector is\n", "($\\frac{1}{\\sqrt{2}}, \\frac{1}{\\sqrt{2}}$) ).\n", "\n", "> Remember that the solution for an eigenvector is not the unique answer;\n", "it only represents a *direction* for an eigenvector corresponding to a\n", "given eigenvalue.\n", "\n", "What are the eigenvectors for the matrix $B$ from the previous section?\n", "\n", "$$B = \\left[\n", " \\begin{array}{ccc} 3 & 2 & 4 \\\\ 2 & 0 & 2 \\\\ 4 & 2 & 3\n", " \\end{array} \\right]$$\n", "\n", "The solution to this problem is available [here](lab3_file/quizzes/eigvec/eigvec.html)\n", "\n", "After solving the questions by hand, you can use Python to check your\n", "answer.\n", "\n", "
\n", "🐾 \n", "
\n", "\n", "[Numpy and Python with Matrices](#Numpy-and-Python-with-Matrices)\n", "\n", "Although the method used here to find the eigenvalues is a direct way to\n", "find the solution, it is not very efficient, especially for large\n", "matrices. Typically, iterative methods such as the Power Method or the\n", "QR algorithm are used (see a linear algebra text such as [Strang (1988)](#Ref:Strang88)\n", "for more details)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Iterative Methods\n", "[lab3:sec:iter]:<#Interative-Methods>\n", "\n", "So far, the only method we’ve seen for solving systems of linear\n", "equations is Gaussian Elimination (with its pivoting variants), which is\n", "only one of a class of *direct methods*. This name derives from the fact\n", "that the Gaussian Elimination algorithm computes the exact solution\n", "*directly*, in a finite number of steps. Other commonly-used direct\n", "methods are based on matrix decomposition or factorizations different\n", "from the $LU$ decomposition (see Section [Decomposition](#Decomposition)); for\n", "example, the $LDL^T$ and Choleski factorizations of a matrix. When the\n", "system to be solved is not too large, it is usually most efficient to\n", "employ a direct technique that minimizes the effects of round-off error\n", "(for example, Gaussian elimination with full pivoting).\n", "\n", "However, the matrices that occur in the discretization of differential\n", "equations are typically *very large* and *sparse* – that is, a large\n", "proportion of the entries in the matrix are zero. In this case, a direct\n", "algorithm, which has a cost on the order of $N^3$ multiplicative\n", "operations, will be spending much of its time inefficiently, operating\n", "on zero entries. In fact, there is another class of solution algorithms\n", "called *iterative methods* which exploit the sparsity of such systems to\n", "reduce the cost of the solution procedure, down to order $N^2$ for\n", "*Jacobi’s method*, the simplest of the iterative methods (see Lab \\#8 )\n", "and as low as order $N$ (the optimal order) for *multigrid methods*\n", "(which we will not discuss here).\n", "\n", "Iterative methods are based on the principle of computing an approximate\n", "solution to a system of equations, where an iterative procedure is\n", "applied to improve the approximation at every iteration. While the exact\n", "answer is never reached, it is hoped that the iterative method will\n", "approach the answer more rapidly than a direct method. For problems\n", "arising from differential equations, this is often possible since these\n", "methods can take advantage of the presence of a large number of zeroes\n", "in the matrix. Even more importantly, most differential equations are\n", "only approximate models of real physical systems in the first place, and\n", "so in many cases, an approximation of the solution is sufficient!!!\n", "\n", "None of the details of iterative methods will be discussed in this Lab.\n", "For now it is enough to know that they exist, and what type of problems\n", "they are used for. Neither will we address the questions: *How quickly\n", "does an iterative method converge to the exact solution?*, *Does it\n", "converge at all?*, and *When are they more efficient than a direct\n", "method?* Iterative methods will be discussed in more detail in Lab \\#8 ,\n", "when a large, sparse system appears in the discretization of a PDE\n", "describing the flow of water in the oceans.\n", "\n", "For even more details on iterative methods, you can also look at [Strang (1988)](#Ref:Strang88) [p. 403ff.], or one of the several books listed in the\n", "Readings section from Lab \\#8 ." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution of an ODE Using Linear Algebra\n", "[lab3:sec:prob]:<#Solution-of-an-ODE-Using-Linear-Algebra>\n", "\n", "So far, we’ve been dealing mainly with matrices with dimensions\n", "$4\\times 4$ at the greatest. If this was the largest linear system that\n", "we ever had to solve, then there would be no need for computers –\n", "everything could be done by hand! Nevertheless, even very simple\n", "differential equations lead to large systems of equations.\n", "\n", "Consider the problem of finding the steady state heat distribution in a\n", "one-dimensional rod, lying along the $x$-axis between 0 and 1. We saw in\n", "Lab \\#1 that the temperature, $u(x)$, can be described by a boundary\n", "value problem, consisting of the ordinary differential equation\n", "$$u_{xx} = f(x),$$ along with boundary values $$u(0)=u(1) = 0.$$ The\n", "only difference between this example and the one from Lab \\#1 is that\n", "the right hand side function, $f(x)$, is non-zero, which corresponds to\n", "a heat source being applied along the length of the rod. The boundary\n", "conditions correspond to the ends of the rod being held at constant\n", "(zero) temperature – this type of condition is known as a fixed or\n", "*Dirichlet* boundary condition.\n", "\n", "If we discretize this equation at $N$ discrete points, $x_i=id$,\n", "$i=0,1,\\dots,N$, where $d = 1/N$ is the grid spacing, then the ordinary\n", "differential equation can be approximated at a point $x_i$ by the\n", "following system of linear equations:\n", "
\n", "(Discrete Differential Equation)\n", "$$\\frac{u_{i+1} - 2u_i+u_{i-1}}{d^2} = f_i,$$ \n", "
\n", "where $f_i=f(x_i)$, and $u_i\\approx u(x_i)$\n", "is an approximation of the steady state temperature at the discrete\n", "points. If we write out all of the equations, for the unknown values\n", "$i=1,\\dots,N-1$, along with the boundary conditions at $i=0,N$, we\n", "obtain the following set of $N+1$ equations in $N+1$ unknowns:\n", "
\n", "(Differential System)\n", "$$\\begin{array}{ccccccccccc}\n", " u_0 & & & & & & & & &=& 0 \\\\\n", " u_0 &-&2u_1 &+& u_2 & & & & &=& d^2 f_1\\\\\n", " & & u_1 &-& 2u_2 &+& u_3 & & &=& d^2f_2\\\\\n", " & & & & & & \\dots & & &=& \\\\\n", " & & & &u_{N-2}&-& 2u_{N-1}&+& u_N &=& d^2f_{N-1}\\\\\n", " & & & & & & & & u_N &=& 0\n", "\\end{array}$$\n", "
\n", "\n", "Remember that this system, like any other linear system, can be written\n", "in matrix notation as\n", "\n", "
\n", "(Differential System Matrix)\n", "$$\\underbrace{\\left[\n", " \\begin{array}{ccccccccc}\n", " 1& 0 & & \\dots & & & & & 0 \\\\\n", " 1& {-2} & {1} & {0} & {\\dots} & && & \\\\\n", " 0& {1} & {-2} & {1} & {0} & {\\dots} & & & \\\\\n", " & {0} & {1} & {-2} & {1} & {0} & {\\dots} & & \\\\\n", " & & & & & & & & \\\\\n", " \\vdots & & & {\\ddots} & {\\ddots} & {\\ddots} & {\\ddots} & {\\ddots} & \\vdots \\\\\n", " & & & & & & & & \\\\\n", " & & & {\\dots} & {0} & {1} & {-2} & {1} & 0 \\\\\n", " & & & &{\\dots} & {0} & {1} & {-2} & 1 \\\\\n", " 0& & & & & \\dots & & 0 & 1 \n", " \\end{array}\n", " \\right]\n", " }_{A_1}\n", " \\underbrace{\\left[\n", " \\begin{array}{c}\n", " u_0 \\\\ {u_1} \\\\ {u_2} \\\\ {u_3} \\\\ \\ \\\\ {\\vdots} \\\\ \\\n", " \\\\ {u_{N-2}} \\\\ {u_{N-1}} \\\\ u_N\n", " \\end{array}\n", " \\right]\n", " }_{U}\n", " = \n", " \\underbrace{\\left[\n", " \\begin{array}{c}\n", " 0 \\\\ {d^2 f_1} \\\\ {d^2 f_2} \\\\ {d^2 f_3} \\\\ \\ \\\\\n", " {\\vdots} \\\\ \\ \\\\ {d^2 f_{N-2}} \\\\ {d^2 f_{N-1}} \\\\ 0 \n", " \\end{array}\n", " \\right] \n", " }_{F}$$\n", "
\n", "\n", "or, simply $A_1 U = F$.\n", "\n", "One question we might ask is: *How well-conditioned is the matrix\n", "$A_1$?* or, in other words, *How easy is this system to solve?* To\n", "answer this question, we leave the right hand side, and consider only\n", "the matrix and its condition number. The size of the condition number is\n", "a measure of how expensive it will be to invert the matrix and hence\n", "solve the discrete boundary value problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem Two\n", "[lab3:prob:dirichlet]:<#Problem-Two> \n", "\n", "> a) Using Python, compute the condition number for\n", "the matrix $A_1$ from Equation [Differential System Matrix](#lab3:eq:dir-system) for several values of $N$\n", "between 5 and 50. ( **Hint:** This will be much easier if you write a\n", "small Python function that outputs the matrix $A$ for a given value of\n", "$N$.)\n", "\n", "> b\\) Can you conjecture how the condition number of $A_1$ depends on $N$?\n", "\n", "> c\\) Another way to write the system of equations is to substitute the\n", "boundary conditions into the equations, and thereby reduce size of the\n", "problem to one of $N-1$ equations in $N-1$ unknowns. The corresponding\n", "matrix is simply the $N-1$ by $N-1$ submatrix of $A_1$\n", "from Equation [Differential System Matrix](#lab3:eq:dir-system) $$A_2 = \\left[\n", " \\begin{array}{ccccccc}\n", " -2 & 1 & 0 & \\dots & && 0 \\\\\n", " 1 & -2 & 1 & 0 & \\dots & & \\\\\n", " 0 & 1 & -2 & 1 & 0 & \\dots & \\\\\n", " & & & & & & \\\\\n", " \\vdots & & \\ddots & \\ddots& \\ddots & \\ddots & \\vdots\\\\\n", " & & & & & & 0 \\\\\n", " & & \\dots & 0 & 1 & -2 & 1 \\\\\n", " 0& & &\\dots & 0 & 1 & -2 \\\\\n", " \\end{array}\n", " \\right]\n", "$$ Does this change in the matrix make a significant difference in the\n", "condition number?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far, we’ve only considered zero Dirichlet boundary values,\n", "$u_0=0=u_N$. Let’s look at a few more types of boundary values …\n", "\n", " **Fixed (non-zero) boundary conditions:**\n", "\n", "> If we fixed the solution at the boundary to be some non-zero values,\n", " say by holding one end at temperature $u_0=a$, and the other at\n", " temperature $u_N=b$, then the matrix itself is not affected. The\n", " only thing that changes in Equation [Differential System](#lab3:eq:dir-system) \n", " is that\n", " a term $a$ is subtracted from the right hand side of the second\n", " equation, and a term $b$ is subtracted from the RHS of the\n", " second-to-last equation. It is clear from what we’ve just said that\n", " non-zero Dirichlet boundary conditions have no effect at all on the\n", " matrix $A_1$ (or $A_2$) since they modify only the right hand side.\n", "\n", " **No-flow boundary conditions:**\n", "\n", "> These are conditions on the first derivative of the temperature\n", " $$u_x(0) = 0,$$ $$u_x(1) = 0,$$ which are also known as *Neumann*\n", " boundary conditions. The requirement that the first derivative of\n", " the temperature be zero at the ends corresponds physically to the\n", " situation where the ends of the rod are *insulated*; that is, rather\n", " than fixing the temperature at the ends of the rod (as we did with\n", " the Dirichlet problem), we require instead that there is no heat\n", " flow in or out of the rod through the ends.\n", "\n", "> There is still one thing that is missing in the mathematical\n", " formulation of this problem: since only derivatives of $u$ appear in\n", " the equations and boundary conditions, the solution is determined\n", " only up to a constant, and for there to be a unique solution, we\n", " must add an extra condition. For example, we could set\n", " $$u \\left(\\frac{1}{2} \\right) = constant,$$ or, more realistically,\n", " say that the total heat contained in the rod is constant, or\n", " $$\\int_0^1 u(x) dx = constant.$$\n", " \n", "> Now, let us look at the discrete formulation of the above problem …\n", "\n", "> The discrete equations do not change, except for that discrete\n", " equations at $i=0,N$ replace the Dirichlet conditions in Equation [Differential System](#lab3:eq:diff-system):\n", "
\n", " (Neumann Boundary Conditions)\n", " $$u_{-1} - 2u_0 +u_{1} = d^2f_0 \\quad {\\rm and} \\quad\n", " u_{N-1} - 2u_N +u_{N+1} = d^2f_N $$ \n", "
\n", " where we have introduced the\n", " additional *ghost points*, or *fictitious points* $u_{-1}$ and\n", " $u_{N+1}$, *lying outside the boundary*. The temperature at these\n", " ghost points can be determined in terms of values in the interior\n", " using the discrete version of the Neumann boundary conditions\n", " $$\\frac{u_{-1} - u_1}{2d} = 0 \\;\\; \\Longrightarrow \\;\\; u_{-1} = u_1,$$\n", " $$\\frac{u_{N+1} - u_{N-1}}{2d} = 0 \\;\\; \\Longrightarrow \\;\\; u_{N+1} = u_{N-1}.$$\n", " Substitute these back into the [Newmann Boundary Conditions](#lab3:eq:neumann-over) to obtain\n", " $$- 2u_0 + 2 u_1 =d^2 f_0 \\quad {\\rm and} \\quad\n", " + 2u_{N-1} - 2 u_N =d^2 f_N .$$ In this case, the matrix is an\n", " $N+1$ by $N+1$ matrix, almost identical to Equation [Differential System Matrix](#lab3:eq:dir-system),\n", " but with the first and last rows slightly modified $$A_3 = \\left[\n", " \\begin{array}{ccccccc}\n", " -2 & 2 & 0 & \\dots & && 0 \\\\\n", " 1 & -2 & 1 & 0 & \\dots & & 0\\\\\n", " 0 & 1 & -2 & 1 & 0 & \\dots & 0\\\\\n", " & & & & & & \\\\\n", " & & & \\ddots& \\ddots & \\ddots & \\\\ \n", " & & & & & & \\\\\n", " 0 & & \\dots & 0 & 1 & -2 & 1 \\\\\n", " 0 & & &\\dots & 0 & 2 & -2\n", " \\end{array}\n", " \\right]$$ \n", " This system is *not solvable*; that is, the $A_3$ above is\n", " singular ( *try it in Python to check for yourself … this should be\n", " easy by modifying the code from [Problem 2](#Problem-Two)).\n", " This is a discrete analogue of the fact that the continuous solution\n", " is not unique. The only way to overcome this problem is to add\n", " another equation for the unknown temperatures." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Physically the reason the problem is not unique is that we don’t know\n", "how hot the rod is. If we think of the full time dependent problem:\n", "\n", "1\\) given fixed temperatures at the end points of the rod (Dirichlet),\n", "whatever the starting temperature of the rod, eventually the rod will\n", "reach equilibrium. with a temperature smoothly varying between the\n", "values (given) at the end points.\n", "\n", "2\\) However, if the rod is insulated (Neumann), no heat can escape and\n", "the final temperature will be related to the initial temperature. To\n", "solve this problem we need to know the steady state,\n", "\n", "a\\) the initial temperature of the rod,\n", "\n", "b\\) the total heat of the rod,\n", "\n", "c\\) or a final temperature at some point of the rod.\n", "\n", "#### Problem Three\n", "\n", "[lab3:prob:neumann]:<#Problem-Three>\n", "\n", "> How can we make the discrete Neumann problem\n", "solvable? Think in terms of discretizing the *solvability conditions*\n", "$u(\\frac{1}{2}) = c$ (condition c) above), or $\\int_0^1 u(x) dx = c$\n", "(condition b) above), (the integral condition can be thought of as an\n", "*average* over the domain, in which case we can approximate it by the\n", "discrete average $\\frac{1}{N}(u_0+u_1+\\dots+u_N)=c$). \n", "\n", "> a) Derive the\n", "matrix corresponding to the linear system to be solved in both of these\n", "cases.\n", "\n", "> b\\) How does the conditioning of the resulting matrix depend on the the\n", "size of the system?\n", "\n", "> c\\) Is it better or worse than for Dirichlet boundary conditions?\n", "\n", " **Periodic boundary conditions:**\n", "\n", "> This refers to the requirement that the temperature at both ends\n", " remains the same: $$u(0) = u(1).$$ Physically, you can think of this\n", " as joining the ends of the rod together, so that it is like a\n", " *ring*. From what we’ve seen already with the other boundary\n", " conditions, it is not hard to see that the discrete form of the\n", " one-dimensional diffusion problem with periodic boundary conditions\n", " leads to an $N\\times N$ matrix of the form $$A_4 = \\left[\n", " \\begin{array}{ccccccc}\n", " -2 & 1 & 0 & \\dots & && 1 \\\\\n", " 1 & -2 & 1 & 0 & \\dots & & 0\\\\\n", " 0 & 1 & -2 & 1 & 0 & \\dots & 0\\\\\n", " & & & & & & \\\\\n", " & & & \\ddots& \\ddots & \\ddots & \\\\ \n", " & & & & & & \\\\\n", " 0 & & \\dots & 0 & 1 & -2 & 1 \\\\\n", " 1 & & &\\dots & 0 & 1 & -2\n", " \\end{array}\n", " \\right],\n", " $$ where the unknown temperatures are now $u_i$, $i=0,1,\\dots, N-1$.\n", " The major change to the form of the matrix is that the elements in\n", " the upper right and lower left corners are now 1 instead of 0. Again\n", " the same problem of the invertibility of the matrix comes up. This\n", " is a symptom of the fact that the continuous problem does not have a\n", " unique solution. It can also be remedied by tacking on an extra\n", " condition, such as in the Neumann problem above.\n", "\n", "#### Problem Four\n", "--------------\n", "[lab3:prob:periodic]: <#Problem-Four> \n", "\n", "> a) Derive the matrix $A_4$ above using the discrete\n", "form [Discrete Differential Equation](#lab3:eq:diff-ode) of the differential equation and the periodic\n", "boundary condition.\n", "\n", "> b) For the periodic problem (with the extra integral condition on the\n", "temperature) how does the conditioning of the matrix compare to that for\n", "the other two discrete problems?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Summary\n", "\n", "As you will have found in these problems, the boundary conditions can\n", "have an influence on the conditioning of a discrete problem.\n", "Furthermore, the method of discretizing the boundary conditions may or\n", "may not have a large effect on the condition number. Consequently, we\n", "must take care when discretizing a problem in order to obtain an\n", "efficient numerical scheme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "
\n", "Strang, G., 1986: Introduction to Applied Mathematics. Wellesley-Cambridge Press, Wellesley, MA.\n", "
\n", "
\n", "Strang, G., 1988: Linear Algebra and its Applications. Harcourt Brace Jovanovich, San Diego, CA, 2nd edition.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numpy and Python with Matrices\n", "\n", "**To start, import numpy,**\n", "\n", "Enter:\n", "\n", "> import numpy as np\n", "\n", "**To enter a matrix,** \n", "$$ A = \\left[ \\begin{array}{ccc} a, & b, & c \\\\\n", " d, & e, & f \\end{array} \\right] $$\n", "Enter:\n", "\n", "> A = np.array([[a, b, c], [d, e, f]])\n", "\n", "**To add two matrices,**\n", "$$ C = A + B$$\n", "\n", "Enter:\n", "\n", "> C = A + B\n", "\n", "**To multiply two matrices,**\n", "$$ C = A \\cdot B $$\n", "\n", "Enter:\n", "\n", "> C = np.dot(A, B)\n", "\n", "**To find the tranpose of a matrix,**\n", "$$ C = A^{T} $$\n", "\n", "Enter:\n", "\n", "> C = A.tranpose()\n", "\n", "**To find the condition number of a matrix,**\n", "\n", "> K = np.linalg.cond(A)\n", "\n", "**To find the inverse of a matrix,**\n", "$$ C = A^{-1} $$\n", "\n", "Enter:\n", "\n", "> C = np.linalg.inv(A)\n", "\n", "**To find the determinant of a matrix,**\n", "$$ K = |A|$$\n", "\n", "Enter:\n", "\n", "> K = np.linalg.det(A)\n", "\n", "**To find the eigenvalues of a matrix,**\n", "\n", "Enter:\n", "\n", "> lamb = np.linalg.eigvals(A)\n", "\n", "**To find the eigenvalues (lamb) and eigenvectors (x) of a matrix,**\n", "\n", "Enter:\n", "\n", "> lamb, x = np.linalg.eig(A)\n", "\n", "lamb[i] are the eigenvalues and x[:, i] are the eigenvectors.\n", "\n", "**To print a matrix,**\n", "$$C$$\n", "\n", "Enter:\n", "\n", "> print (C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Glossary\n", "[glossary.unnumbered]:<#Glossary>\n", "\n", "**A** \n", "\n", "augmented matrix\n", "\n", "> The $m \\times (n+1)$ matrix representing a linear system,\n", " $Ax = b$, with the right hand side vector appended to the\n", " coefficient matrix: $$\\left[ \n", " \\begin{array}{cc}\n", " \\begin{array}{ccccc} \n", " a_{11} & & \\ldots & & a_{1n} \\\\\n", " \\vdots & & \\ddots & & \\vdots \\\\\n", " a_{m1} & & \\ldots & & a_{mn} \n", " \\end{array}\n", " &\n", " \\left| \n", " \\begin{array}{rc}\n", " & b_{1} \\\\ & \\vdots \\\\ & b_{m}\n", " \\end{array} \n", " \\right. \n", " \\end{array}\n", " \\right]$$\n", "\n", "> The right most column is the right hand side vector or augmented\n", " column.\n", "\n", "**C** \n", "\n", "characteristic equation\n", "\n", "> The equation:\n", " $$\\det(A - \\lambda I) = 0 , \\ \\ \\ \\ or \\ \\ \\ \\ Ax = \\lambda x$$\n", " where $A$ is a *square matrix*, $I$ is the *identity matrix*,\n", " $\\lambda$ is an *eigenvalue* of $A$, and $x$ is the corresponding\n", " *eigenvector* of $\\lambda$.\n", "\n", "coefficient matrix\n", "\n", "> A $m \\times n$ matrix made up with the coefficients $a_{ij}$ of the\n", " $n$ unknowns from the $m$ equations of a set of linear equations,\n", " where $i$ is the row index and $j$ is the column index: $$\\left[\n", " \\begin{array}{ccccccc}\n", " & a_{11} & & \\ldots & & a_{1n} & \\\\\n", " & \\vdots & & \\ddots & & \\vdots & \\\\\n", " & a_{m1} & & \\ldots & & a_{mn} &\n", " \\end{array}\n", " \\right]$$\n", "\n", "condition number\n", "\n", "> A number, $K$, that refers to the sensitivity of a *nonsingular*\n", " matrix, $A$, i.e. given a system $Ax = b$, $K$ reflects whether\n", " small changes in $A$ and $b$ will have any effect on the solution.\n", " The matrix is well-conditioned if $K$ is close to one. The number is\n", " described as: $$K(A) = \\|A\\| \\|A^{-1}\\| \n", " \\ \\ \\ \\ or \\ \\ \\ \\\n", " K(A) = \\frac{\\lambda_{max}}{\\lambda_{min}}$$ where $\\lambda_{max}$\n", " and $\\lambda_{min}$ are largest and smallest *eigenvalues* of $A$\n", " respectively.\n", "\n", "**D** \n", "\n", "decomposition\n", "\n", "> Factoring a matrix, $A$, into two factors, e.g., the Gaussian\n", " elimination amounts to factoring $A$ into a product of two matrices.\n", " One is the lower triangular matrix, $L$, and the other is the upper\n", " triangular matrix, $U$.\n", "\n", "diagonal matrix\n", "\n", "> A square matrix with the entries $a_{ij} = 0 $ whenever $i \\neq j$.\n", "\n", "**E**\n", "\n", "eigenvalue\n", "\n", "> A number, $\\lambda$, that must satisfy the *characteristic equation*\n", " $\\det(A - \\lambda I) = 0.$\n", "\n", "eigenvector\n", "\n", "> A vector, $x$, which corresponds to an *eigenvalue* of a *square\n", " matrix* $A$, satisfying the characteristic equation:\n", " $$Ax = \\lambda x .$$\n", "\n", "**H** \n", "\n", "homogeneous equations\n", "\n", "> A set of linear equations, $Ax = b$ with the zero vector on the\n", " right hand side, i.e. $b = 0$.\n", "\n", "**I** \n", "\n", "inhomogeneous equations\n", "\n", "> A set of linear equations, $Ax = b$ such that $b \\neq 0$.\n", "\n", "identity matrix\n", "\n", "> A *diagonal matrix* with the entries $a_{ii} = 1$:\n", " $$\\left[ \\begin{array}{ccccccc}\n", " & 1 & 0 & \\ldots & \\ldots & 0 & \\\\\n", " & 0 & 1 & \\ddots & & \\vdots & \\\\\n", " & \\vdots & \\ddots & \\ddots & \\ddots & \\vdots \\\\\n", " & \\vdots & & \\ddots & 1 & 0 & \\\\\n", " & 0 & \\ldots & \\ldots & 0 & 1 &\n", " \\end{array} \\right]$$\n", "\n", "ill-conditioned matrix\n", "\n", "> A matrix with a large *condition number*, i.e., the matrix is not\n", " well-behaved, and small errors to the matrix will have great effects\n", " to the solution.\n", "\n", "invertible matrix\n", "\n", "> A square matrix, $A$, such that there exists another matrix,\n", " $A^{-1}$, which satisfies:\n", " $$AA^{-1} = I \\ \\ \\ \\ and \\ \\ \\ \\ A^{-1}A = I$$\n", "\n", "> The matrix, $A^{-1}$, is the *inverse* of $A$. An invertible matrix\n", " is *nonsingular*.\n", "\n", "**L** \n", "\n", "linear system\n", "\n", "> A set of $m$ equations in $n$ unknowns: $$\\begin{array}{ccccccc}\n", " a_{11}x_{1} & + & \\ldots & + & a_{1n}x_{n} & = & b_{1} \\\\\n", " a_{21}x_{1} & + & \\ldots & + & a_{2n}x_{n} & = & b_{2} \\\\\n", " & & \\vdots & & & & \\vdots \\\\\n", " a_{m1}x_{1} & + & \\ldots & + & a_{mn}x_{n} & = & b_{m} \n", " \\end{array}$$ with unknowns $x_{i}$ and coefficients\n", " $a_{ij}, b_{j}$.\n", "\n", "lower triangular matrix\n", "\n", "> A square matrix, $L$, with the entries $l_{ij} = 0$, whenever\n", " $j > i$: $$\\left[\n", " \\begin{array}{ccccccc}\n", " & * & 0 & \\ldots & \\ldots & 0 & \\\\\n", " & * & * & \\ddots & & \\vdots & \\\\\n", " & \\vdots & & \\ddots & \\ddots & \\vdots & \\\\\n", " & \\vdots & & & * & 0 & \\\\\n", " & * & \\ldots & \\ldots & \\ldots & * &\n", " \\end{array}\n", " \\right]$$\n", "\n", "**N** \n", "\n", "nonsingular matrix\n", "\n", "> A square matrix,$A$, that is invertible, i.e. the system $Ax = b$\n", " has a *unique solution*.\n", "\n", "**S** \n", "\n", "singular matrix\n", "\n", "> A $n \\times n$ matrix that is degenerate and does not have an\n", " inverse (refer to *invertible*), i.e., the system $Ax = b$ does not\n", " have a *unique solution*.\n", "\n", "sparse matrix\n", "\n", "> A matrix with a high percentage of zero entries.\n", "\n", "square matrix\n", "\n", "> A matrix with the same number of rows and columns.\n", "\n", "**T** \n", "\n", "transpose\n", "\n", "> A $n \\times m$ matrix, $A^{T}$, that has the columns of a\n", " $m \\times n$ matrix, $A$, as its rows, and the rows of $A$ as its\n", " columns, i.e. the entry $a_{ij}$ in $A$ becomes $a_{ji}$ in $A^{T}$,\n", " e.g.\n", "\n", "> $$A = \n", " \\left[ \\begin{array}{ccc} 1 & 2 & 3 \\\\ 4 & 5 & 6 \\end{array} \\right] \n", " \\ \\ \\rightarrow \\ \\ \n", " A^{T} = \n", " \\left[ \\begin{array}{cc} 1 & 4 \\\\ 2 & 5 \\\\ 3 & 6 \\end{array} \\right]$$\n", "\n", "tridiagonal matrix\n", "\n", "> A square matrix with the entries $a_{ij} = 0$, $| i-j | > 1 $:\n", "\n", "> $$\\left[\n", " \\begin{array}{cccccccc}\n", " & * & * & 0 & \\ldots & \\ldots & 0 & \\\\\n", " & * & * & * & \\ddots & & \\vdots & \\\\\n", " & 0 & * & \\ddots & \\ddots & \\ddots & \\vdots & \\\\\n", " & \\vdots & \\ddots & \\ddots & \\ddots & * & 0 & \\\\\n", " & \\vdots & & \\ddots & * & * & * & \\\\\n", " & 0 & \\ldots & \\ldots & 0 & * & * &\n", " \\end{array}\n", " \\right]$$\n", "\n", "**U** \n", "\n", "unique solution\n", "\n", "> There is only solution, $x$, that satisfies a particular linear\n", " system, $Ax = b$, for the given $A$. That is, this linear system has\n", " exactly one solution. The matrix $A$ of the system is *invertible*\n", " or *nonsingular*.\n", "\n", "upper triangular matrix\n", "\n", "> A square matrix, $U$, with the entries $u_{ij} = 0$ whenever\n", " $i > j$:\n", "\n", "> $$\\left[\n", " \\begin{array}{ccccccc}\n", " & * & \\ldots & \\ldots & \\ldots & * & \\\\\n", " & 0 & * & & & \\vdots & \\\\\n", " & \\vdots & \\ddots & \\ddots & & \\vdots & \\\\\n", " & \\vdots & & \\ddots & * & * & \\\\\n", " & 0 & \\ldots & \\ldots & 0 & * &\n", " \\end{array}\n", " \\right]$$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "header_metadata": { "chead": "Feb., 2020", "lhead": "Numeric Lab 2" }, "jupytext": { "cell_metadata_filter": "all", "encoding": "# -*- coding: utf-8 -*-", "formats": "ipynb,py:percent", "notebook_metadata_filter": "all,-language_info,-toc,-latex_envs" }, "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.6" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "meta-9" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "324.176px" }, "toc_section_display": "block", "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }