diff --git a/README.md b/README.md index 73ca14a..c7011f5 100644 --- a/README.md +++ b/README.md @@ -79,10 +79,10 @@ general, I will not post homework solutions. | |2/16|8|Linear Algebra| |6|2/21|9|Linear systems: Gauss elimination| | |2/23|10|Linear Systems: LU factorization| -|7|2/28||Midterm Review| -| |3/2||Midterm| -|8|3/7|11|Linear Systems: Error analysis| -| |3/9|13|Eigenvalues| +|7|2/28|11|Linear Systems: Error analysis| +| |3/2|12|Eigenvalues| +|8|3/7|1-10 |Midterm Review| +| |3/9|1-10|Midterm| |9|3/14| N/A |Spring Break!| | |3/16| N/A |Spring Break!| |10|3/21|12|Linear Systems: Iterative methods| diff --git a/lecture_09/mass_springs.png b/lecture_09/mass_springs.png new file mode 100644 index 0000000..d0dc65c Binary files /dev/null and b/lecture_09/mass_springs.png differ diff --git a/lecture_11/.ipynb_checkpoints/lecture_11-checkpoint.ipynb b/lecture_11/.ipynb_checkpoints/lecture_11-checkpoint.ipynb new file mode 100644 index 0000000..e2aa425 --- /dev/null +++ b/lecture_11/.ipynb_checkpoints/lecture_11-checkpoint.ipynb @@ -0,0 +1,699 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What is the determinant of A?\n", + "\n", + "$A=\\left[ \\begin{array}{ccc}\n", + "10 & 2 & 1 \\\\\n", + "0 & 1 & 1 \\\\\n", + "0 & 0 & 10\\end{array} \\right]$\n", + "\n", + "![responses to determinant of A](det_A.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# LU Decomposition\n", + "### efficient storage of matrices for solutions\n", + "\n", + "Considering the same solution set:\n", + "\n", + "$y=Ax$\n", + "\n", + "Assume that we can perform Gauss elimination and achieve this formula:\n", + "\n", + "$Ux=d$ \n", + "\n", + "Where, $U$ is an upper triangular matrix that we derived from Gauss elimination and $d$ is the set of dependent variables after Gauss elimination. \n", + "\n", + "Assume there is a lower triangular matrix, $L$, with ones on the diagonal and same dimensions of $U$ and the following is true:\n", + "\n", + "$L(Ux-d)=Ax-y=0$\n", + "\n", + "Now, $Ax=LUx$, so $A=LU$, and $y=Ld$.\n", + "\n", + "$2x_{1}+x_{2}=1$\n", + "\n", + "$x_{1}+3x_{2}=1$\n", + "\n", + "\n", + "$\\left[ \\begin{array}{cc}\n", + "2 & 1 \\\\\n", + "1 & 3 \\end{array} \\right]\n", + "\\left[\\begin{array}{c} \n", + "x_{1} \\\\ \n", + "x_{2} \\end{array}\\right]=\n", + "\\left[\\begin{array}{c} \n", + "1 \\\\\n", + "1\\end{array}\\right]$\n", + "\n", + "f21=0.5\n", + "\n", + "A(2,1)=1-1 = 0 \n", + "\n", + "A(2,2)=3-0.5=2.5\n", + "\n", + "y(2)=1-0.5=0.5\n", + "\n", + "$L(Ux-d)=\n", + "\\left[ \\begin{array}{cc}\n", + "1 & 0 \\\\\n", + "0.5 & 1 \\end{array} \\right]\n", + "\\left(\\left[ \\begin{array}{cc}\n", + "2 & 1 \\\\\n", + "0 & 2.5 \\end{array} \\right]\n", + "\\left[\\begin{array}{c} \n", + "x_{1} \\\\ \n", + "x_{2} \\end{array}\\right]-\n", + "\\left[\\begin{array}{c} \n", + "1 \\\\\n", + "0.5\\end{array}\\right]\\right)=0$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000\n", + " 0.50000 1.00000\n", + "\n", + "U =\n", + "\n", + " 2.00000 1.00000\n", + " 0.00000 2.50000\n", + "\n", + "ans =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "d =\n", + "\n", + " 1.00000\n", + " 0.50000\n", + "\n", + "y =\n", + "\n", + " 1\n", + " 1\n", + "\n" + ] + } + ], + "source": [ + "A=[2,1;1,3]\n", + "L=[1,0;0.5,1]\n", + "U=[2,1;0,2.5]\n", + "L*U\n", + "\n", + "d=[1;0.5]\n", + "y=L*d" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pivoting for LU factorization\n", + "\n", + "LU factorization uses the same method as Gauss elimination so it is also necessary to perform partial pivoting when creating the lower and upper triangular matrices. \n", + "\n", + "Matlab and Octave use pivoting in the command \n", + "\n", + "`[L,U,P]=lu(A)`\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n", + "\n", + " -- Built-in Function: [L, U] = lu (A)\n", + " -- Built-in Function: [L, U, P] = lu (A)\n", + " -- Built-in Function: [L, U, P, Q] = lu (S)\n", + " -- Built-in Function: [L, U, P, Q, R] = lu (S)\n", + " -- Built-in Function: [...] = lu (S, THRES)\n", + " -- Built-in Function: Y = lu (...)\n", + " -- Built-in Function: [...] = lu (..., \"vector\")\n", + " Compute the LU decomposition of A.\n", + "\n", + " If A is full subroutines from LAPACK are used and if A is sparse\n", + " then UMFPACK is used.\n", + "\n", + " The result is returned in a permuted form, according to the\n", + " optional return value P. For example, given the matrix 'a = [1, 2;\n", + " 3, 4]',\n", + "\n", + " [l, u, p] = lu (A)\n", + "\n", + " returns\n", + "\n", + " l =\n", + "\n", + " 1.00000 0.00000\n", + " 0.33333 1.00000\n", + "\n", + " u =\n", + "\n", + " 3.00000 4.00000\n", + " 0.00000 0.66667\n", + "\n", + " p =\n", + "\n", + " 0 1\n", + " 1 0\n", + "\n", + " The matrix is not required to be square.\n", + "\n", + " When called with two or three output arguments and a spare input\n", + " matrix, 'lu' does not attempt to perform sparsity preserving column\n", + " permutations. Called with a fourth output argument, the sparsity\n", + " preserving column transformation Q is returned, such that 'P * A *\n", + " Q = L * U'.\n", + "\n", + " Called with a fifth output argument and a sparse input matrix, 'lu'\n", + " attempts to use a scaling factor R on the input matrix such that 'P\n", + " * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n", + " stable factorization.\n", + "\n", + " An additional input argument THRES, that defines the pivoting\n", + " threshold can be given. THRES can be a scalar, in which case it\n", + " defines the UMFPACK pivoting tolerance for both symmetric and\n", + " unsymmetric cases. If THRES is a 2-element vector, then the first\n", + " element defines the pivoting tolerance for the unsymmetric UMFPACK\n", + " pivoting strategy and the second for the symmetric strategy. By\n", + " default, the values defined by 'spparms' are used ([0.1, 0.001]).\n", + "\n", + " Given the string argument \"vector\", 'lu' returns the values of P\n", + " and Q as vector values, such that for full matrix, 'A (P,:) = L *\n", + " U', and 'R(P,:) * A (:, Q) = L * U'.\n", + "\n", + " With two output arguments, returns the permuted forms of the upper\n", + " and lower triangular matrices, such that 'A = L * U'. With one\n", + " output argument Y, then the matrix returned by the LAPACK routines\n", + " is returned. If the input matrix is sparse then the matrix L is\n", + " embedded into U to give a return value similar to the full case.\n", + " For both full and sparse matrices, 'lu' loses the permutation\n", + " information.\n", + "\n", + " See also: luupdate, ilu, chol, hess, qr, qz, schur, svd.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help lu" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5e-05\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.00015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.00025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0003\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tLU decomp\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tLU decomp\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tOctave \\\\\n", + "\n", + "\t\n", + "\t\tOctave \\\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% time LU solution vs backslash\n", + "t_lu=zeros(100,1);\n", + "t_bs=zeros(100,1);\n", + "for N=1:100\n", + " A=rand(N,N);\n", + " y=rand(N,1);\n", + " [L,U,P]=lu(A);\n", + "\n", + " tic; d=L\\y; x=U\\d; t_lu(N)=toc;\n", + "\n", + " tic; x=A\\y; t_bs(N)=toc;\n", + "end\n", + "plot([1:100],t_lu,[1:100],t_bs) \n", + "legend('LU decomp','Octave \\\\')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n", + "\n", + "![Springs-masses](../lecture_09/mass_springs.svg)\n", + "\n", + "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n", + "\n", + "$m_{1}g+k(x_{2}-x_{1})-kx_{1}=0$\n", + "\n", + "$m_{2}g+k(x_{3}-x_{2})-k(x_{2}-x_{1})=0$\n", + "\n", + "$m_{3}g+k(x_{4}-x_{3})-k(x_{3}-x_{2})=0$\n", + "\n", + "$m_{4}g-k(x_{4}-x_{3})=0$\n", + "\n", + "in matrix form:\n", + "\n", + "$\\left[ \\begin{array}{cccc}\n", + "2k & -k & 0 & 0 \\\\\n", + "-k & 2k & -k & 0 \\\\\n", + "0 & -k & 2k & -k \\\\\n", + "0 & 0 & -k & k \\end{array} \\right]\n", + "\\left[ \\begin{array}{c}\n", + "x_{1} \\\\\n", + "x_{2} \\\\\n", + "x_{3} \\\\\n", + "x_{4} \\end{array} \\right]=\n", + "\\left[ \\begin{array}{c}\n", + "m_{1}g \\\\\n", + "m_{2}g \\\\\n", + "m_{3}g \\\\\n", + "m_{4}g \\end{array} \\right]$" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 20 -10 0 0\n", + " -10 20 -10 0\n", + " 0 -10 20 -10\n", + " 0 0 -10 10\n", + "\n", + "y =\n", + "\n", + " 9.8100\n", + " 19.6200\n", + " 29.4300\n", + " 39.2400\n", + "\n" + ] + } + ], + "source": [ + "k=10; % N/m\n", + "m1=1; % kg\n", + "m2=2;\n", + "m3=3;\n", + "m4=4;\n", + "g=9.81; % m/s^2\n", + "K=[2*k -k 0 0; -k 2*k -k 0; 0 -k 2*k -k; 0 0 -k k]\n", + "y=[m1*g;m2*g;m3*g;m4*g]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This matrix, K, is symmetric. \n", + "\n", + "`K(i,j)==K(j,i)`\n", + "\n", + "Now we can use,\n", + "\n", + "## Cholesky Factorization\n", + "\n", + "We can decompose the matrix, K into two matrices, $U$ and $U^{T}$, where\n", + "\n", + "$K=U^{T}U$\n", + "\n", + "each of the components of U can be calculated with the following equations:\n", + "\n", + "$u_{ii}=\\sqrt{a_{ii}-\\sum_{k=1}^{i-1}u_{ki}^{2}}$\n", + "\n", + "$u_{ij}=\\frac{a_{ij}-\\sum_{k=1}^{i-1}u_{ki}u_{kj}}{u_{ii}}$\n", + "\n", + "so for K" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 20 -10 0 0\n", + " -10 20 -10 0\n", + " 0 -10 20 -10\n", + " 0 0 -10 10\n", + "\n" + ] + } + ], + "source": [ + "K" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "u11 = 4.4721\n", + "u12 = -2.2361\n", + "u13 = 0\n", + "u14 = 0\n", + "u22 = 3.8730\n", + "u23 = -2.5820\n", + "u24 = 0\n", + "u33 = 3.6515\n", + "u34 = -2.7386\n", + "u44 = 1.5811\n", + "U =\n", + "\n", + " 4.47214 -2.23607 0.00000 0.00000\n", + " 0.00000 3.87298 -2.58199 0.00000\n", + " 0.00000 0.00000 3.65148 -2.73861\n", + " 0.00000 0.00000 0.00000 1.58114\n", + "\n" + ] + } + ], + "source": [ + "u11=sqrt(K(1,1))\n", + "u12=(K(1,2))/u11\n", + "u13=(K(1,3))/u11\n", + "u14=(K(1,4))/u11\n", + "u22=sqrt(K(2,2)-u12^2)\n", + "u23=(K(2,3)-u12*u13)/u22\n", + "u24=(K(2,4)-u12*u14)/u22\n", + "u33=sqrt(K(3,3)-u13^2-u23^2)\n", + "u34=(K(3,4)-u13*u14-u23*u24)/u33\n", + "u44=sqrt(K(4,4)-u14^2-u24^2-u34^2)\n", + "U=[u11,u12,u13,u14;0,u22,u23,u24;0,0,u33,u34;0,0,0,u44]" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans =\n", + "\n", + " 20.00000 -10.00000 0.00000 0.00000\n", + " -10.00000 20.00000 -10.00000 0.00000\n", + " 0.00000 -10.00000 20.00000 -10.00000\n", + " 0.00000 0.00000 -10.00000 10.00000\n", + "\n" + ] + } + ], + "source": [ + "U'*U" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "average time spent for Cholesky factored solution = 1.623154e-05+/-1.166726e-05\n", + "average time spent for backslash solution = 1.675844e-05+/-1.187234e-05\n" + ] + } + ], + "source": [ + "% time solution for Cholesky vs backslash\n", + "t_chol=zeros(1000,1);\n", + "t_bs=zeros(1000,1);\n", + "for i=1:1000\n", + " tic; d=U'*y; x=U\\d; t_chol(i)=toc;\n", + " tic; x=K\\y; t_bs(i)=toc;\n", + "end\n", + "fprintf('average time spent for Cholesky factored solution = %e+/-%e',mean(t_chol),std(t_chol))\n", + "\n", + "fprintf('average time spent for backslash solution = %e+/-%e',mean(t_bs),std(t_bs))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_11/det_A.png b/lecture_11/det_A.png new file mode 100644 index 0000000..50f6ac1 Binary files /dev/null and b/lecture_11/det_A.png differ diff --git a/lecture_11/lecture_11.aux b/lecture_11/lecture_11.aux index cf3c530..a9c4771 100644 --- a/lecture_11/lecture_11.aux +++ b/lecture_11/lecture_11.aux @@ -17,12 +17,22 @@ \providecommand\HyField@AuxAddToFields[1]{} \providecommand\HyField@AuxAddToCoFields[2]{} \providecommand \oddpage@label [2]{} -\@writefile{toc}{\contentsline {section}{\numberline {1}LU Decomposition}{1}{section.1}} -\newlabel{lu-decomposition}{{1}{1}{LU Decomposition}{section.1}{}} -\@writefile{toc}{\contentsline {subsubsection}{\numberline {1.0.1}efficient storage of matrices for solutions}{1}{subsubsection.1.0.1}} -\newlabel{efficient-storage-of-matrices-for-solutions}{{1.0.1}{1}{efficient storage of matrices for solutions}{subsubsection.1.0.1}{}} -\@writefile{toc}{\contentsline {subsection}{\numberline {1.1}Pivoting for LU factorization}{2}{subsection.1.1}} -\newlabel{pivoting-for-lu-factorization}{{1.1}{2}{Pivoting for LU factorization}{subsection.1.1}{}} -\@writefile{lof}{\contentsline {figure}{\numberline {1}{\ignorespaces Springs-masses\relax }}{5}{figure.caption.1}} -\@writefile{toc}{\contentsline {subsection}{\numberline {1.2}Cholesky Factorization}{6}{subsection.1.2}} -\newlabel{cholesky-factorization}{{1.2}{6}{Cholesky Factorization}{subsection.1.2}{}} +\@writefile{toc}{\contentsline {subsection}{\numberline {0.1}My question from last class}{1}{subsection.0.1}} +\newlabel{my-question-from-last-class}{{0.1}{1}{My question from last class}{subsection.0.1}{}} +\@writefile{lof}{\contentsline {figure}{\numberline {1}{\ignorespaces responses to determinant of A\relax }}{1}{figure.caption.1}} +\@writefile{toc}{\contentsline {subsection}{\numberline {0.2}Your questions from last class}{1}{subsection.0.2}} +\newlabel{your-questions-from-last-class}{{0.2}{1}{Your questions from last class}{subsection.0.2}{}} +\@writefile{toc}{\contentsline {subsection}{\numberline {0.3}Midterm preference}{2}{subsection.0.3}} +\newlabel{midterm-preference}{{0.3}{2}{Midterm preference}{subsection.0.3}{}} +\@writefile{lof}{\contentsline {figure}{\numberline {2}{\ignorespaces responses to midterm date\relax }}{2}{figure.caption.2}} +\@writefile{toc}{\contentsline {subsubsection}{\numberline {0.3.1}Midterm Questions}{2}{subsubsection.0.3.1}} +\newlabel{midterm-questions}{{0.3.1}{2}{Midterm Questions}{subsubsection.0.3.1}{}} +\@writefile{toc}{\contentsline {section}{\numberline {1}LU Decomposition}{2}{section.1}} +\newlabel{lu-decomposition}{{1}{2}{LU Decomposition}{section.1}{}} +\@writefile{toc}{\contentsline {subsubsection}{\numberline {1.0.1}efficient storage of matrices for solutions}{2}{subsubsection.1.0.1}} +\newlabel{efficient-storage-of-matrices-for-solutions}{{1.0.1}{2}{efficient storage of matrices for solutions}{subsubsection.1.0.1}{}} +\@writefile{toc}{\contentsline {subsection}{\numberline {1.1}Pivoting for LU factorization}{4}{subsection.1.1}} +\newlabel{pivoting-for-lu-factorization}{{1.1}{4}{Pivoting for LU factorization}{subsection.1.1}{}} +\@writefile{lof}{\contentsline {figure}{\numberline {3}{\ignorespaces Springs-masses\relax }}{7}{figure.caption.3}} +\@writefile{toc}{\contentsline {subsection}{\numberline {1.2}Cholesky Factorization}{8}{subsection.1.2}} +\newlabel{cholesky-factorization}{{1.2}{8}{Cholesky Factorization}{subsection.1.2}{}} diff --git a/lecture_11/lecture_11.ipynb b/lecture_11/lecture_11.ipynb index c290a47..41f1870 100644 --- a/lecture_11/lecture_11.ipynb +++ b/lecture_11/lecture_11.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 18, + "execution_count": 1, "metadata": { "collapsed": true }, @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 2, "metadata": { "collapsed": true }, @@ -22,6 +22,51 @@ "setdefaults" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## My question from last class \n", + "\n", + "$A=\\left[ \\begin{array}{ccc}\n", + "10 & 2 & 1 \\\\\n", + "0 & 1 & 1 \\\\\n", + "0 & 0 & 10\\end{array} \\right]$\n", + "\n", + "![responses to determinant of A](det_A.png)\n", + "\n", + "## Your questions from last class\n", + "\n", + "1. Need more linear algebra review\n", + " \n", + " -We will keep doing Linear Algebra, try the practice problems in 'linear_algebra'\n", + "\n", + "2. How do I do HW3? \n", + " \n", + " -demo today\n", + "\n", + "3. For hw4 is the spring constant (K) suppose to be given? \n", + " \n", + " -yes, its 30 N/m\n", + " \n", + "4. Deapool or Joker?\n", + "\n", + "\n", + "## Midterm preference\n", + "\n", + "![responses to midterm date](midterm_date.png)\n", + "\n", + "### Midterm Questions\n", + "\n", + "1. Notes allowed\n", + " \n", + " -no\n", + "\n", + "2. Will there be a review/study sheet\n", + "\n", + " -yes" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -137,6 +182,35 @@ "y=L*d" ] }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 5.0000\n", + "ans = 1\n", + "ans = 5\n", + "ans = 5\n", + "ans = 5.0000\n" + ] + } + ], + "source": [ + "% what is the determinant of L, U and A?\n", + "\n", + "det(A)\n", + "det(L)\n", + "det(U)\n", + "det(L)*det(U)\n", + "det(L*U)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -250,7 +324,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -297,68 +371,58 @@ "\n", "\n", "\t\n", - "\t\t\n", + "\t\t\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t5e-05\n", + "\t\t\n", + "\t\t0.0005\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0.0001\n", + "\t\t\n", + "\t\t0.001\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0.00015\n", + "\t\t\n", + "\t\t0.0015\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0.0002\n", + "\t\t\n", + "\t\t0.002\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0.00025\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.0003\n", - "\t\n", - "\n", - "\n", - "\t\t\n", + "\t\t\n", "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t20\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t40\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t60\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t80\n", "\t\n", "\n", @@ -370,7 +434,7 @@ "\n", "\n", "\n", - "\t\n", + "\t\n", "\n", "\n", "\n", @@ -386,7 +450,7 @@ "\t\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", "\tOctave \\\\\n", "\n", @@ -395,7 +459,7 @@ "\t\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -425,9 +489,9 @@ " y=rand(N,1);\n", " [L,U,P]=lu(A);\n", "\n", - " tic; d=L\\y; x=U\\d; t_lu(N)=toc;\n", + " tic; d=inv(L)*y; x=inv(U)*d; t_lu(N)=toc;\n", "\n", - " tic; x=A\\y; t_bs(N)=toc;\n", + " tic; x=inv(A)*y; t_bs(N)=toc;\n", "end\n", "plot([1:100],t_lu,[1:100],t_bs) \n", "legend('LU decomp','Octave \\\\')" @@ -474,7 +538,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 10, "metadata": { "collapsed": false }, @@ -538,7 +602,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -563,7 +627,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 12, "metadata": { "collapsed": false }, @@ -608,7 +672,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -619,21 +683,21 @@ "text": [ "ans =\n", "\n", - " 20.00000 -10.00000 0.00000 0.00000\n", - " -10.00000 20.00000 -10.00000 0.00000\n", - " 0.00000 -10.00000 20.00000 -10.00000\n", - " 0.00000 0.00000 -10.00000 10.00000\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + " 1 1 1 1\n", "\n" ] } ], "source": [ - "U'*U" + "(U'*U)'==U'*U" ] }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -642,8 +706,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "average time spent for Cholesky factored solution = 1.623154e-05+/-1.166726e-05\n", - "average time spent for backslash solution = 1.675844e-05+/-1.187234e-05\n" + "average time spent for Cholesky factored solution = 1.465964e-05+/-9.806001e-06\n", + "average time spent for backslash solution = 1.555967e-05+/-1.048561e-05\n" ] } ], @@ -659,6 +723,15 @@ "\n", "fprintf('average time spent for backslash solution = %e+/-%e',mean(t_bs),std(t_bs))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/lecture_11/lecture_11.log b/lecture_11/lecture_11.log index 10423a7..981cd40 100644 --- a/lecture_11/lecture_11.log +++ b/lecture_11/lecture_11.log @@ -1,4 +1,4 @@ -This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=pdflatex 2017.1.11) 22 FEB 2017 22:01 +This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=pdflatex 2017.1.11) 27 FEB 2017 11:55 entering extended mode restricted \write18 enabled. %&-line parsing enabled. @@ -734,76 +734,89 @@ LaTeX Font Info: Try loading font information for T1+cmtt on input line 279. (/usr/share/texlive/texmf-dist/tex/latex/base/t1cmtt.fd File: t1cmtt.fd 2014/09/29 v2.5h Standard LaTeX font definitions ) -LaTeX Font Info: Font shape `T1/ppl/bx/n' in size <14.4> not available -(Font) Font shape `T1/ppl/b/n' tried instead on input line 287. -LaTeX Font Info: Font shape `T1/ppl/bx/n' in size <10.95> not available -(Font) Font shape `T1/ppl/b/n' tried instead on input line 290. - - -Package hyperref Warning: Difference (2) between bookmark levels is greater -(hyperref) than one, level fixed on input line 290. - +LaTeX Font Info: Font shape `T1/ppl/bx/n' in size <12> not available +(Font) Font shape `T1/ppl/b/n' tried instead on input line 288. LaTeX Font Info: Font shape `U/msa/m/n' will be -(Font) scaled to size 11.40997pt on input line 294. +(Font) scaled to size 11.40997pt on input line 290. LaTeX Font Info: Font shape `U/msa/m/n' will be -(Font) scaled to size 8.33606pt on input line 294. +(Font) scaled to size 8.33606pt on input line 290. LaTeX Font Info: Font shape `U/msa/m/n' will be -(Font) scaled to size 6.25204pt on input line 294. +(Font) scaled to size 6.25204pt on input line 290. LaTeX Font Info: Font shape `U/msb/m/n' will be -(Font) scaled to size 11.40997pt on input line 294. +(Font) scaled to size 11.40997pt on input line 290. LaTeX Font Info: Font shape `U/msb/m/n' will be -(Font) scaled to size 8.33606pt on input line 294. +(Font) scaled to size 8.33606pt on input line 290. LaTeX Font Info: Font shape `U/msb/m/n' will be -(Font) scaled to size 6.25204pt on input line 294. -[1 +(Font) scaled to size 6.25204pt on input line 290. + + +File: det_A.png Graphic file (type png) + +Package pdftex.def Info: det_A.png used on input line 294. +(pdftex.def) Requested size: 375.80544pt x 202.3582pt. + [1 + +{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map} <./det_A.png>] + +File: midterm_date.png Graphic file (type png) + +Package pdftex.def Info: midterm_date.png used on input line 324. +(pdftex.def) Requested size: 375.80544pt x 171.71808pt. +LaTeX Font Info: Font shape `T1/ppl/bx/n' in size <10.95> not available +(Font) Font shape `T1/ppl/b/n' tried instead on input line 328. +LaTeX Font Info: Font shape `T1/ppl/bx/n' in size <14.4> not available +(Font) Font shape `T1/ppl/b/n' tried instead on input line 342. -{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map}] -LaTeX Font Info: Font shape `T1/ppl/bx/n' in size <12> not available -(Font) Font shape `T1/ppl/b/n' tried instead on input line 372. -LaTeX Font Info: Try loading font information for TS1+cmtt on input line 387 + +Package hyperref Warning: Difference (2) between bookmark levels is greater +(hyperref) than one, level fixed on input line 345. + +[2 <./midterm_date.png>] [3] +LaTeX Font Info: Try loading font information for TS1+cmtt on input line 461 . (/usr/share/texlive/texmf-dist/tex/latex/base/ts1cmtt.fd File: ts1cmtt.fd 2014/09/29 v2.5h Standard LaTeX font definitions -) [2] [3] +) [4] [5] LaTeX Font Info: Font shape `T1/cmtt/bx/n' in size <10.95> not available -(Font) Font shape `T1/cmtt/m/n' tried instead on input line 473. +(Font) Font shape `T1/cmtt/m/n' tried instead on input line 547. - -File: lecture_11_files/lecture_11_6_0.pdf Graphic file (type pdf) + +File: lecture_11_files/lecture_11_8_0.pdf Graphic file (type pdf) - -Package pdftex.def Info: lecture_11_files/lecture_11_6_0.pdf used on input line - 487. + +Package pdftex.def Info: lecture_11_files/lecture_11_8_0.pdf used on input line + 561. (pdftex.def) Requested size: 449.6789pt x 337.25917pt. - [4] -Underfull \hbox (badness 10000) in paragraph at lines 489--490 + +Underfull \hbox (badness 10000) in paragraph at lines 563--564 [] - -File: mass_springs.png Graphic file (type png) - -Package pdftex.def Info: mass_springs.png used on input line 497. -(pdftex.def) Requested size: 89.93611pt x 152.576pt. +<../lecture_09/mass_springs.png, id=80, 112.42pt x 190.7125pt> +File: ../lecture_09/mass_springs.png Graphic file (type png) -[5 <./lecture_11_files/lecture_11_6_0.pdf> <./mass_springs.png>] [6] [7] -Package atveryend Info: Empty hook `BeforeClearDocument' on input line 656. + +Package pdftex.def Info: ../lecture_09/mass_springs.png used on input line 571. -[8] -Package atveryend Info: Empty hook `AfterLastShipout' on input line 656. +(pdftex.def) Requested size: 89.93611pt x 152.576pt. + [6 <./lecture_11_files/lecture_11_8_0.pdf>] [7 <../lecture_09/mass_springs.png +>] [8] +Package atveryend Info: Empty hook `BeforeClearDocument' on input line 734. + [9] +Package atveryend Info: Empty hook `AfterLastShipout' on input line 734. (./lecture_11.aux) -Package atveryend Info: Executing hook `AtVeryEndDocument' on input line 656. -Package atveryend Info: Executing hook `AtEndAfterFileList' on input line 656. +Package atveryend Info: Executing hook `AtVeryEndDocument' on input line 734. +Package atveryend Info: Executing hook `AtEndAfterFileList' on input line 734. Package rerunfilecheck Info: File `lecture_11.out' has not changed. -(rerunfilecheck) Checksum: 06B2288B56AA590C3216AAE67971C54B;297. -Package atveryend Info: Empty hook `AtVeryVeryEnd' on input line 656. +(rerunfilecheck) Checksum: 41743ED7C91140C56E53FBD51DB19D65;568. +Package atveryend Info: Empty hook `AtVeryVeryEnd' on input line 734. ) Here is how much of TeX's memory you used: - 10920 strings out of 493029 - 163179 string characters out of 6136234 - 268780 words of memory out of 5000000 - 14206 multiletter control sequences out of 15000+600000 + 10955 strings out of 493029 + 163773 string characters out of 6136234 + 268999 words of memory out of 5000000 + 14228 multiletter control sequences out of 15000+600000 39052 words of font info for 103 fonts, out of 8000000 for 9000 1141 hyphenation exceptions out of 8191 36i,9n,77p,460b,465s stack positions out of 5000i,500n,10000p,200000b,80000s @@ -817,10 +830,10 @@ hare/texmf/fonts/type1/public/cm-super/sfit1095.pfb> -Output written on lecture_11.pdf (8 pages, 154124 bytes). +Output written on lecture_11.pdf (9 pages, 191934 bytes). PDF statistics: - 117 PDF objects out of 1000 (max. 8388607) - 92 compressed objects within 1 object stream - 14 named destinations out of 1000 (max. 500000) - 43 words of extra memory for PDF output out of 10000 (max. 10000000) + 152 PDF objects out of 1000 (max. 8388607) + 121 compressed objects within 2 object streams + 27 named destinations out of 1000 (max. 500000) + 85 words of extra memory for PDF output out of 10000 (max. 10000000) diff --git a/lecture_11/lecture_11.out b/lecture_11/lecture_11.out index 41c2cd3..481e738 100644 --- a/lecture_11/lecture_11.out +++ b/lecture_11/lecture_11.out @@ -1,4 +1,8 @@ -\BOOKMARK [1][-]{section.1}{LU Decomposition}{}% 1 -\BOOKMARK [2][-]{subsubsection.1.0.1}{efficient storage of matrices for solutions}{section.1}% 2 -\BOOKMARK [2][-]{subsection.1.1}{Pivoting for LU factorization}{section.1}% 3 -\BOOKMARK [2][-]{subsection.1.2}{Cholesky Factorization}{section.1}% 4 +\BOOKMARK [2][-]{subsection.0.1}{My question from last class}{}% 1 +\BOOKMARK [2][-]{subsection.0.2}{Your questions from last class}{}% 2 +\BOOKMARK [2][-]{subsection.0.3}{Midterm preference}{}% 3 +\BOOKMARK [3][-]{subsubsection.0.3.1}{Midterm Questions}{subsection.0.3}% 4 +\BOOKMARK [1][-]{section.1}{LU Decomposition}{}% 5 +\BOOKMARK [2][-]{subsubsection.1.0.1}{efficient storage of matrices for solutions}{section.1}% 6 +\BOOKMARK [2][-]{subsection.1.1}{Pivoting for LU factorization}{section.1}% 7 +\BOOKMARK [2][-]{subsection.1.2}{Cholesky Factorization}{section.1}% 8 diff --git a/lecture_11/lecture_11.pdf b/lecture_11/lecture_11.pdf index 3e2c828..9fa1342 100644 Binary files a/lecture_11/lecture_11.pdf and b/lecture_11/lecture_11.pdf differ diff --git a/lecture_11/lecture_11.tex b/lecture_11/lecture_11.tex index 22d4064..318852d 100644 --- a/lecture_11/lecture_11.tex +++ b/lecture_11/lecture_11.tex @@ -277,13 +277,68 @@ \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}18}]:} \PY{c}{\PYZpc{}plot \PYZhy{}\PYZhy{}format svg} +{\color{incolor}In [{\color{incolor}1}]:} \PY{c}{\PYZpc{}plot \PYZhy{}\PYZhy{}format svg} \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}19}]:} \PY{n}{setdefaults} +{\color{incolor}In [{\color{incolor}2}]:} \PY{n}{setdefaults} \end{Verbatim} + \subsection{My question from last +class}\label{my-question-from-last-class} + +\(A=\left[ \begin{array}{ccc} 10 & 2 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 10\end{array} \right]\) + +\begin{figure}[htbp] +\centering +\includegraphics{det_A.png} +\caption{responses to determinant of A} +\end{figure} + +\subsection{Your questions from last +class}\label{your-questions-from-last-class} + +\begin{enumerate} +\def\labelenumi{\arabic{enumi}.} +\item + Need more linear algebra review + + -We will keep doing Linear Algebra, try the practice problems in + 'linear\_algebra' +\item + How do I do HW3? + + -demo today +\item + For hw4 is the spring constant (K) suppose to be given? + + -yes, its 30 N/m +\item + Deapool or Joker? +\end{enumerate} + +\subsection{Midterm preference}\label{midterm-preference} + +\begin{figure}[htbp] +\centering +\includegraphics{midterm_date.png} +\caption{responses to midterm date} +\end{figure} + +\subsubsection{Midterm Questions}\label{midterm-questions} + +\begin{enumerate} +\def\labelenumi{\arabic{enumi}.} +\item + Notes allowed + + -no +\item + Will there be a review/study sheet + + -yes +\end{enumerate} + \section{LU Decomposition}\label{lu-decomposition} \subsubsection{efficient storage of matrices for @@ -366,6 +421,25 @@ \subsubsection{efficient storage of matrices for 1 + \end{Verbatim} + + \begin{Verbatim}[commandchars=\\\{\}] +{\color{incolor}In [{\color{incolor}5}]:} \PY{c}{\PYZpc{} what is the determinant of L, U and A?} + + \PY{n+nb}{det}\PY{p}{(}\PY{n}{A}\PY{p}{)} + \PY{n+nb}{det}\PY{p}{(}\PY{n}{L}\PY{p}{)} + \PY{n+nb}{det}\PY{p}{(}\PY{n}{U}\PY{p}{)} + \PY{n+nb}{det}\PY{p}{(}\PY{n}{L}\PY{p}{)}\PY{o}{*}\PY{n+nb}{det}\PY{p}{(}\PY{n}{U}\PY{p}{)} + \PY{n+nb}{det}\PY{p}{(}\PY{n}{L}\PY{o}{*}\PY{n}{U}\PY{p}{)} +\end{Verbatim} + + \begin{Verbatim}[commandchars=\\\{\}] +ans = 5.0000 +ans = 1 +ans = 5 +ans = 5 +ans = 5.0000 + \end{Verbatim} \subsection{Pivoting for LU @@ -467,24 +541,24 @@ \subsubsection{efficient storage of matrices for \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}22}]:} \PY{c}{\PYZpc{} time LU solution vs backslash} - \PY{n}{t\PYZus{}lu}\PY{p}{=}\PY{n+nb}{zeros}\PY{p}{(}\PY{l+m+mi}{100}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} - \PY{n}{t\PYZus{}bs}\PY{p}{=}\PY{n+nb}{zeros}\PY{p}{(}\PY{l+m+mi}{100}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} - \PY{k}{for} \PY{n}{N}\PY{p}{=}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{100} - \PY{n}{A}\PY{p}{=}\PY{n+nb}{rand}\PY{p}{(}\PY{n}{N}\PY{p}{,}\PY{n}{N}\PY{p}{)}\PY{p}{;} - \PY{n}{y}\PY{p}{=}\PY{n+nb}{rand}\PY{p}{(}\PY{n}{N}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} - \PY{p}{[}\PY{n}{L}\PY{p}{,}\PY{n}{U}\PY{p}{,}\PY{n}{P}\PY{p}{]}\PY{p}{=}\PY{n+nb}{lu}\PY{p}{(}\PY{n}{A}\PY{p}{)}\PY{p}{;} - - \PY{n+nb}{tic}\PY{p}{;} \PY{n}{d}\PY{p}{=}\PY{n}{L}\PY{o}{\PYZbs{}}\PY{n}{y}\PY{p}{;} \PY{n}{x}\PY{p}{=}\PY{n}{U}\PY{o}{\PYZbs{}}\PY{n}{d}\PY{p}{;} \PY{n}{t\PYZus{}lu}\PY{p}{(}\PY{n}{N}\PY{p}{)}\PY{p}{=}\PY{n+nb}{toc}\PY{p}{;} - - \PY{n+nb}{tic}\PY{p}{;} \PY{n}{x}\PY{p}{=}\PY{n}{A}\PY{o}{\PYZbs{}}\PY{n}{y}\PY{p}{;} \PY{n}{t\PYZus{}bs}\PY{p}{(}\PY{n}{N}\PY{p}{)}\PY{p}{=}\PY{n+nb}{toc}\PY{p}{;} - \PY{k}{end} - \PY{n+nb}{plot}\PY{p}{(}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{100}\PY{p}{]}\PY{p}{,}\PY{n}{t\PYZus{}lu}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{100}\PY{p}{]}\PY{p}{,}\PY{n}{t\PYZus{}bs}\PY{p}{)} - \PY{n+nb}{legend}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{LU decomp\PYZsq{}}\PY{p}{,}\PY{l+s}{\PYZsq{}}\PY{l+s}{Octave \PYZbs{}\PYZbs{}\PYZsq{}}\PY{p}{)} +{\color{incolor}In [{\color{incolor}9}]:} \PY{c}{\PYZpc{} time LU solution vs backslash} + \PY{n}{t\PYZus{}lu}\PY{p}{=}\PY{n+nb}{zeros}\PY{p}{(}\PY{l+m+mi}{100}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} + \PY{n}{t\PYZus{}bs}\PY{p}{=}\PY{n+nb}{zeros}\PY{p}{(}\PY{l+m+mi}{100}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} + \PY{k}{for} \PY{n}{N}\PY{p}{=}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{100} + \PY{n}{A}\PY{p}{=}\PY{n+nb}{rand}\PY{p}{(}\PY{n}{N}\PY{p}{,}\PY{n}{N}\PY{p}{)}\PY{p}{;} + \PY{n}{y}\PY{p}{=}\PY{n+nb}{rand}\PY{p}{(}\PY{n}{N}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} + \PY{p}{[}\PY{n}{L}\PY{p}{,}\PY{n}{U}\PY{p}{,}\PY{n}{P}\PY{p}{]}\PY{p}{=}\PY{n+nb}{lu}\PY{p}{(}\PY{n}{A}\PY{p}{)}\PY{p}{;} + + \PY{n+nb}{tic}\PY{p}{;} \PY{n}{d}\PY{p}{=}\PY{n+nb}{inv}\PY{p}{(}\PY{n}{L}\PY{p}{)}\PY{o}{*}\PY{n}{y}\PY{p}{;} \PY{n}{x}\PY{p}{=}\PY{n+nb}{inv}\PY{p}{(}\PY{n}{U}\PY{p}{)}\PY{o}{*}\PY{n}{d}\PY{p}{;} \PY{n}{t\PYZus{}lu}\PY{p}{(}\PY{n}{N}\PY{p}{)}\PY{p}{=}\PY{n+nb}{toc}\PY{p}{;} + + \PY{n+nb}{tic}\PY{p}{;} \PY{n}{x}\PY{p}{=}\PY{n+nb}{inv}\PY{p}{(}\PY{n}{A}\PY{p}{)}\PY{o}{*}\PY{n}{y}\PY{p}{;} \PY{n}{t\PYZus{}bs}\PY{p}{(}\PY{n}{N}\PY{p}{)}\PY{p}{=}\PY{n+nb}{toc}\PY{p}{;} + \PY{k}{end} + \PY{n+nb}{plot}\PY{p}{(}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{100}\PY{p}{]}\PY{p}{,}\PY{n}{t\PYZus{}lu}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{100}\PY{p}{]}\PY{p}{,}\PY{n}{t\PYZus{}bs}\PY{p}{)} + \PY{n+nb}{legend}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{LU decomp\PYZsq{}}\PY{p}{,}\PY{l+s}{\PYZsq{}}\PY{l+s}{Octave \PYZbs{}\PYZbs{}\PYZsq{}}\PY{p}{)} \end{Verbatim} \begin{center} - \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{lecture_11_files/lecture_11_6_0.pdf} + \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{lecture_11_files/lecture_11_8_0.pdf} \end{center} { \hspace*{\fill} \\} @@ -494,7 +568,7 @@ \subsubsection{efficient storage of matrices for \begin{figure}[htbp] \centering -\includegraphics{mass_springs.png} +\includegraphics{../lecture_09/mass_springs.png} \caption{Springs-masses} \end{figure} @@ -514,7 +588,7 @@ \subsubsection{efficient storage of matrices for \(\left[ \begin{array}{cccc} 2k & -k & 0 & 0 \\ -k & 2k & -k & 0 \\ 0 & -k & 2k & -k \\ 0 & 0 & -k & k \end{array} \right] \left[ \begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{array} \right]= \left[ \begin{array}{c} m_{1}g \\ m_{2}g \\ m_{3}g \\ m_{4}g \end{array} \right]\) \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}24}]:} \PY{n}{k}\PY{p}{=}\PY{l+m+mi}{10}\PY{p}{;} \PY{c}{\PYZpc{} N/m} +{\color{incolor}In [{\color{incolor}10}]:} \PY{n}{k}\PY{p}{=}\PY{l+m+mi}{10}\PY{p}{;} \PY{c}{\PYZpc{} N/m} \PY{n}{m1}\PY{p}{=}\PY{l+m+mi}{1}\PY{p}{;} \PY{c}{\PYZpc{} kg} \PY{n}{m2}\PY{p}{=}\PY{l+m+mi}{2}\PY{p}{;} \PY{n}{m3}\PY{p}{=}\PY{l+m+mi}{3}\PY{p}{;} @@ -565,7 +639,7 @@ \subsection{Cholesky Factorization}\label{cholesky-factorization} so for K \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}25}]:} \PY{n}{K} +{\color{incolor}In [{\color{incolor}11}]:} \PY{n}{K} \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] @@ -580,7 +654,7 @@ \subsection{Cholesky Factorization}\label{cholesky-factorization} \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}26}]:} \PY{n}{u11}\PY{p}{=}\PY{n+nb}{sqrt}\PY{p}{(}\PY{n}{K}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{)} +{\color{incolor}In [{\color{incolor}12}]:} \PY{n}{u11}\PY{p}{=}\PY{n+nb}{sqrt}\PY{p}{(}\PY{n}{K}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{)} \PY{n}{u12}\PY{p}{=}\PY{p}{(}\PY{n}{K}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}\PY{o}{/}\PY{n}{u11} \PY{n}{u13}\PY{p}{=}\PY{p}{(}\PY{n}{K}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{3}\PY{p}{)}\PY{p}{)}\PY{o}{/}\PY{n}{u11} \PY{n}{u14}\PY{p}{=}\PY{p}{(}\PY{n}{K}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{)}\PY{p}{)}\PY{o}{/}\PY{n}{u11} @@ -615,22 +689,22 @@ \subsection{Cholesky Factorization}\label{cholesky-factorization} \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}27}]:} \PY{n}{U}\PY{o}{\PYZsq{}}\PY{o}{*}\PY{n}{U} +{\color{incolor}In [{\color{incolor}17}]:} \PY{p}{(}\PY{n}{U}\PY{o}{\PYZsq{}}\PY{o}{*}\PY{n}{U}\PY{p}{)}\PY{o}{\PYZsq{}}\PY{o}{==}\PY{n}{U}\PY{o}{\PYZsq{}}\PY{o}{*}\PY{n}{U} \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] ans = - 20.00000 -10.00000 0.00000 0.00000 - -10.00000 20.00000 -10.00000 0.00000 - 0.00000 -10.00000 20.00000 -10.00000 - 0.00000 0.00000 -10.00000 10.00000 + 1 1 1 1 + 1 1 1 1 + 1 1 1 1 + 1 1 1 1 \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] -{\color{incolor}In [{\color{incolor}37}]:} \PY{c}{\PYZpc{} time solution for Cholesky vs backslash} +{\color{incolor}In [{\color{incolor}18}]:} \PY{c}{\PYZpc{} time solution for Cholesky vs backslash} \PY{n}{t\PYZus{}chol}\PY{p}{=}\PY{n+nb}{zeros}\PY{p}{(}\PY{l+m+mi}{1000}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} \PY{n}{t\PYZus{}bs}\PY{p}{=}\PY{n+nb}{zeros}\PY{p}{(}\PY{l+m+mi}{1000}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{)}\PY{p}{;} \PY{k}{for} \PY{n}{i}\PY{p}{=}\PY{l+m+mi}{1}\PY{p}{:}\PY{l+m+mi}{1000} @@ -643,11 +717,15 @@ \subsection{Cholesky Factorization}\label{cholesky-factorization} \end{Verbatim} \begin{Verbatim}[commandchars=\\\{\}] -average time spent for Cholesky factored solution = 1.623154e-05+/-1.166726e-05 -average time spent for backslash solution = 1.675844e-05+/-1.187234e-05 +average time spent for Cholesky factored solution = 1.465964e-05+/-9.806001e-06 +average time spent for backslash solution = 1.555967e-05+/-1.048561e-05 \end{Verbatim} + \begin{Verbatim}[commandchars=\\\{\}] +{\color{incolor}In [{\color{incolor} }]:} +\end{Verbatim} + % Add a bibliography block to the postdoc diff --git a/lecture_11/lecture_11_files/lecture_11_8_0.pdf b/lecture_11/lecture_11_files/lecture_11_8_0.pdf new file mode 100644 index 0000000..7f562c6 Binary files /dev/null and b/lecture_11/lecture_11_files/lecture_11_8_0.pdf differ diff --git a/lecture_11/lecture_11_files/lecture_11_8_0.svg b/lecture_11/lecture_11_files/lecture_11_8_0.svg new file mode 100644 index 0000000..37a1c3f --- /dev/null +++ b/lecture_11/lecture_11_files/lecture_11_8_0.svg @@ -0,0 +1,139 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0 + + + + + 20 + + + + + 40 + + + + + 60 + + + + + 80 + + + + + 100 + + + + + + + + + + + + + LU decomp + + + + + LU decomp + + + + + + Octave \\ + + + Octave \ + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_11/midterm_date.png b/lecture_11/midterm_date.png new file mode 100644 index 0000000..4dd0663 Binary files /dev/null and b/lecture_11/midterm_date.png differ diff --git a/lecture_11/octave-workspace b/lecture_11/octave-workspace index 2e1bdfa..8c437bb 100644 Binary files a/lecture_11/octave-workspace and b/lecture_11/octave-workspace differ diff --git a/lecture_12/.ipynb_checkpoints/lecture_12-checkpoint.ipynb b/lecture_12/.ipynb_checkpoints/lecture_12-checkpoint.ipynb new file mode 100644 index 0000000..1332b29 --- /dev/null +++ b/lecture_12/.ipynb_checkpoints/lecture_12-checkpoint.ipynb @@ -0,0 +1,1483 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 47, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## My question from last class \n", + "\n", + "\n", + "## Your questions from last class\n", + "\n", + "1. Need more linear algebra review\n", + " \n", + " -We will keep doing Linear Algebra, try the practice problems in 'linear_algebra'\n", + "\n", + "2. How do I do HW3? \n", + " \n", + " -demo today\n", + "\n", + "3. For hw4 is the spring constant (K) suppose to be given? \n", + " \n", + " -yes, its 30 N/m\n", + " \n", + "4. Deapool or Joker?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Matrix Inverse and Condition\n", + "\n", + "\n", + "Considering the same solution set:\n", + "\n", + "$y=Ax$\n", + "\n", + "If we know that $A^{-1}A=I$, then \n", + "\n", + "$A^{-1}y=A^{-1}Ax=x$\n", + "\n", + "so \n", + "\n", + "$x=A^{-1}y$\n", + "\n", + "Where, $A^{-1}$ is the inverse of matrix $A$.\n", + "\n", + "$2x_{1}+x_{2}=1$\n", + "\n", + "$x_{1}+3x_{2}=1$\n", + "\n", + "$Ax=y$\n", + "\n", + "$\\left[ \\begin{array}{cc}\n", + "2 & 1 \\\\\n", + "1 & 3 \\end{array} \\right]\n", + "\\left[\\begin{array}{c} \n", + "x_{1} \\\\ \n", + "x_{2} \\end{array}\\right]=\n", + "\\left[\\begin{array}{c} \n", + "1 \\\\\n", + "1\\end{array}\\right]$\n", + "\n", + "$A^{-1}=\\frac{1}{2*3-1*1}\\left[ \\begin{array}{cc}\n", + "3 & 1 \\\\\n", + "-1 & 2 \\end{array} \\right]=\n", + "\\left[ \\begin{array}{cc}\n", + "3/5 & -1/5 \\\\\n", + "-1/5 & 2/5 \\end{array} \\right]$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "invA =\n", + "\n", + " 0.60000 -0.20000\n", + " -0.20000 0.40000\n", + "\n", + "ans =\n", + "\n", + " 1.00000 0.00000\n", + " 0.00000 1.00000\n", + "\n", + "ans =\n", + "\n", + " 1.00000 0.00000\n", + " 0.00000 1.00000\n", + "\n" + ] + } + ], + "source": [ + "A=[2,1;1,3]\n", + "invA=1/5*[3,-1;-1,2]\n", + "\n", + "A*invA\n", + "invA*A" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "How did we know the inverse of A? \n", + "\n", + "for 2$\\times$2 matrices, it is always:\n", + "\n", + "$A=\\left[ \\begin{array}{cc}\n", + "A_{11} & A_{12} \\\\\n", + "A_{21} & A_{22} \\end{array} \\right]$\n", + "\n", + "$A^{-1}=\\frac{1}{det(A)}\\left[ \\begin{array}{cc}\n", + "A_{22} & -A_{12} \\\\\n", + "-A_{21} & A_{11} \\end{array} \\right]$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$AA^{-1}=\\frac{1}{A_{11}A_{22}-A_{21}A_{12}}\\left[ \\begin{array}{cc}\n", + "A_{11}A_{22}-A_{21}A_{12} & -A_{11}A_{12}+A_{12}A_{11} \\\\\n", + "A_{21}A_{22}-A_{22}A_{21} & -A_{21}A_{12}+A_{22}A_{11} \\end{array} \\right]\n", + "=\\left[ \\begin{array}{cc}\n", + "1 & 0 \\\\\n", + "0 & 1 \\end{array} \\right]$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What about bigger matrices?\n", + "\n", + "We can use the LU-decomposition\n", + "\n", + "$A=LU$\n", + "\n", + "$A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$\n", + "\n", + "if we divide $A^{-1}$ into n-column vectors, $a_{n}$, then\n", + "\n", + "$Aa_{1}=\\left[\\begin{array}{c} \n", + "1 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$Aa_{2}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "1 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$Aa_{n}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "1 \\end{array} \\right]$\n", + "\n", + "\n", + "Which we can solve for each $a_{n}$ with LU-decomposition, knowing the lower and upper triangular decompositions, then \n", + "\n", + "$A^{-1}=\\left[ \\begin{array}{cccc}\n", + "| & | & & | \\\\\n", + "a_{1} & a_{2} & \\cdots & a_{3} \\\\\n", + "| & | & & | \\end{array} \\right]$\n", + "\n", + "\n", + "$Ld_{1}=\\left[\\begin{array}{c} \n", + "1 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$;~Ua_{1}=d_{1}$\n", + "\n", + "$Ld_{2}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "1 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$;~Ua_{2}=d_{2}$\n", + "\n", + "$Ld_{n}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "1 \\\\ \n", + "\\vdots \\\\\n", + "n \\end{array} \\right]$\n", + "$;~Ua_{n}=d_{n}$\n", + "\n", + "Consider the following matrix:\n", + "\n", + "$A=\\left[ \\begin{array}{ccc}\n", + "2 & -1 & 0\\\\\n", + "-1 & 2 & -1\\\\\n", + "0 & -1 & 1 \\end{array} \\right]$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 -1 0\n", + " -1 2 -1\n", + " 0 -1 1\n", + "\n", + "U =\n", + "\n", + " 2.00000 -1.00000 0.00000\n", + " 0.00000 1.50000 -1.00000\n", + " 0.00000 -1.00000 1.00000\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " -0.50000 1.00000 0.00000\n", + " 0.00000 0.00000 1.00000\n", + "\n" + ] + } + ], + "source": [ + "A=[2,-1,0;-1,2,-1;0,-1,1]\n", + "U=A;\n", + "L=eye(3,3);\n", + "U(2,:)=U(2,:)-U(2,1)/U(1,1)*U(1,:)\n", + "L(2,1)=A(2,1)/A(1,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "L =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " -0.50000 1.00000 0.00000\n", + " 0.00000 -0.66667 1.00000\n", + "\n", + "U =\n", + "\n", + " 2.00000 -1.00000 0.00000\n", + " 0.00000 1.50000 -1.00000\n", + " 0.00000 0.00000 0.33333\n", + "\n" + ] + } + ], + "source": [ + "L(3,2)=U(3,2)/U(2,2)\n", + "U(3,:)=U(3,:)-U(3,2)/U(2,2)*U(2,:)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now solve for $d_1$ then $a_1$, $d_2$ then $a_2$, and $d_3$ then $a_{3}$\n", + "\n", + "$Ld_{1}=\\left[\\begin{array}{c} \n", + "1 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$;~Ua_{1}=d_{1}$" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "d1 =\n", + "\n", + " 1.00000\n", + " 0.50000\n", + " 0.33333\n", + "\n" + ] + } + ], + "source": [ + "d1=zeros(3,1);\n", + "d1(1)=1;\n", + "d1(2)=0-L(2,1)*d1(1);\n", + "d1(3)=0-L(3,1)*d1(1)-L(3,2)*d1(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a1 =\n", + "\n", + " 1.00000\n", + " 1.00000\n", + " 1.00000\n", + "\n" + ] + } + ], + "source": [ + "a1=zeros(3,1);\n", + "a1(3)=d1(3)/U(3,3);\n", + "a1(2)=1/U(2,2)*(d1(2)-U(2,3)*a1(3));\n", + "a1(1)=1/U(1,1)*(d1(1)-U(1,2)*a1(2)-U(1,3)*a1(3))" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "d2 =\n", + "\n", + " 0.00000\n", + " 1.00000\n", + " 0.66667\n", + "\n" + ] + } + ], + "source": [ + "d2=zeros(3,1);\n", + "d2(1)=0;\n", + "d2(2)=1-L(2,1)*d2(1);\n", + "d2(3)=0-L(3,1)*d2(1)-L(3,2)*d2(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a2 =\n", + "\n", + " 1.0000\n", + " 2.0000\n", + " 2.0000\n", + "\n" + ] + } + ], + "source": [ + "a2=zeros(3,1);\n", + "a2(3)=d2(3)/U(3,3);\n", + "a2(2)=1/U(2,2)*(d2(2)-U(2,3)*a2(3));\n", + "a2(1)=1/U(1,1)*(d2(1)-U(1,2)*a2(2)-U(1,3)*a2(3))" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "d3 =\n", + "\n", + " 0\n", + " 0\n", + " 1\n", + "\n" + ] + } + ], + "source": [ + "d3=zeros(3,1);\n", + "d3(1)=0;\n", + "d3(2)=0-L(2,1)*d3(1);\n", + "d3(3)=1-L(3,1)*d3(1)-L(3,2)*d3(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a3 =\n", + "\n", + " 1.00000\n", + " 2.00000\n", + " 3.00000\n", + "\n" + ] + } + ], + "source": [ + "a3=zeros(3,1);\n", + "a3(3)=d3(3)/U(3,3);\n", + "a3(2)=1/U(2,2)*(d3(2)-U(2,3)*a3(3));\n", + "a3(1)=1/U(1,1)*(d3(1)-U(1,2)*a3(2)-U(1,3)*a3(3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Final solution for $A^{-1}$ is $[a_{1}~a_{2}~a_{3}]$" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "invA =\n", + "\n", + " 1.00000 1.00000 1.00000\n", + " 1.00000 2.00000 2.00000\n", + " 1.00000 2.00000 3.00000\n", + "\n", + "ans =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " 0.00000 1.00000 -0.00000\n", + " -0.00000 -0.00000 1.00000\n", + "\n" + ] + } + ], + "source": [ + "invA=[a1,a2,a3]\n", + "A*invA" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now the solution of $x$ to $Ax=y$ is $x=A^{-1}y$" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "y =\n", + "\n", + " 1\n", + " 2\n", + " 3\n", + "\n", + "x =\n", + "\n", + " 6.0000\n", + " 11.0000\n", + " 14.0000\n", + "\n", + "xbs =\n", + "\n", + " 6.0000\n", + " 11.0000\n", + " 14.0000\n", + "\n", + "ans =\n", + "\n", + " -3.5527e-15\n", + " -8.8818e-15\n", + " -1.0658e-14\n", + "\n", + "ans = 2.2204e-16\n" + ] + } + ], + "source": [ + "y=[1;2;3]\n", + "x=invA*y\n", + "xbs=A\\y\n", + "x-xbs\n", + "eps" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tinversion\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tinversion\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbackslash\n", + "\n", + "\t\n", + "\t\tbackslash\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmultiplication\n", + "\n", + "\t\n", + "\t\tmultiplication\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=100;\n", + "n=[1:N];\n", + "t_inv=zeros(N,1);\n", + "t_bs=zeros(N,1);\n", + "t_mult=zeros(N,1);\n", + "for i=1:N\n", + " A=rand(i,i);\n", + " tic\n", + " invA=inv(A);\n", + " t_inv(i)=toc;\n", + " b=rand(i,1);\n", + " tic;\n", + " x=A\\b;\n", + " t_bs(i)=toc;\n", + " tic;\n", + " x=invA*b;\n", + " t_mult(i)=toc;\n", + "end\n", + "plot(n,t_inv,n,t_bs,n,t_mult)\n", + "axis([0 100 0 0.002])\n", + "legend('inversion','backslash','multiplication','Location','NorthWest')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Condition of a matrix \n", + "### *just checked in to see what condition my condition was in*\n", + "### Matrix norms\n", + "\n", + "The Euclidean norm of a vector is measure of the magnitude (in 3D this would be: $|x|=\\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$) in general the equation is:\n", + "\n", + "$||x||_{e}=\\sqrt{\\sum_{i=1}^{n}x_{i}^{2}}$\n", + "\n", + "For a matrix, A, the same norm is called the Frobenius norm:\n", + "\n", + "$||A||_{f}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{2}}$\n", + "\n", + "In general we can calculate any $p$-norm where\n", + "\n", + "$||A||_{p}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{p}}$\n", + "\n", + "so the p=1, 1-norm is \n", + "\n", + "$||A||_{1}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{1}}=\\sum_{i=1}^{n}\\sum_{i=1}^{m}|A_{i,j}|$\n", + "\n", + "$||A||_{\\infty}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{\\infty}}=\\max_{1\\le i \\le n}\\sum_{j=1}^{m}|A_{i,j}|$\n", + "\n", + "### Condition of Matrix\n", + "\n", + "The matrix condition is the product of \n", + "\n", + "$Cond(A) = ||A||\\cdot||A^{-1}||$ \n", + "\n", + "So each norm will have a different condition number, but the limit is $Cond(A)\\ge 1$\n", + "\n", + "An estimate of the rounding error is based on the condition of A:\n", + "\n", + "$\\frac{||\\Delta x||}{x} \\le Cond(A) \\frac{||\\Delta A||}{||A||}$\n", + "\n", + "So if the coefficients of A have accuracy to $10^{-t}\n", + "\n", + "and the condition of A, $Cond(A)=10^{c}$\n", + "\n", + "then the solution for x can have rounding errors up to $10^{c-t}$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 1.00000 0.50000 0.33333\n", + " 0.50000 0.33333 0.25000\n", + " 0.33333 0.25000 0.20000\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " 0.50000 1.00000 0.00000\n", + " 0.33333 1.00000 1.00000\n", + "\n", + "U =\n", + "\n", + " 1.00000 0.50000 0.33333\n", + " 0.00000 0.08333 0.08333\n", + " 0.00000 -0.00000 0.00556\n", + "\n" + ] + } + ], + "source": [ + "A=[1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5]\n", + "[L,U]=LU_naive(A)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "invA=" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000\n", + " 0.50000 1.00000\n", + "\n", + "U =\n", + "\n", + " 2.00000 1.00000\n", + " 0.00000 2.50000\n", + "\n", + "ans =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "d =\n", + "\n", + " 1.00000\n", + " 0.50000\n", + "\n", + "y =\n", + "\n", + " 1\n", + " 1\n", + "\n" + ] + } + ], + "source": [ + "A=[2,1;1,3]\n", + "L=[1,0;0.5,1]\n", + "U=[2,1;0,2.5]\n", + "L*U\n", + "\n", + "d=[1;0.5]\n", + "y=L*d" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 5.0000\n", + "ans = 1\n", + "ans = 5\n", + "ans = 5\n", + "ans = 5.0000\n" + ] + } + ], + "source": [ + "% what is the determinant of L, U and A?\n", + "\n", + "det(A)\n", + "det(L)\n", + "det(U)\n", + "det(L)*det(U)\n", + "det(L*U)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pivoting for LU factorization\n", + "\n", + "LU factorization uses the same method as Gauss elimination so it is also necessary to perform partial pivoting when creating the lower and upper triangular matrices. \n", + "\n", + "Matlab and Octave use pivoting in the command \n", + "\n", + "`[L,U,P]=lu(A)`\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n", + "\n", + " -- Built-in Function: [L, U] = lu (A)\n", + " -- Built-in Function: [L, U, P] = lu (A)\n", + " -- Built-in Function: [L, U, P, Q] = lu (S)\n", + " -- Built-in Function: [L, U, P, Q, R] = lu (S)\n", + " -- Built-in Function: [...] = lu (S, THRES)\n", + " -- Built-in Function: Y = lu (...)\n", + " -- Built-in Function: [...] = lu (..., \"vector\")\n", + " Compute the LU decomposition of A.\n", + "\n", + " If A is full subroutines from LAPACK are used and if A is sparse\n", + " then UMFPACK is used.\n", + "\n", + " The result is returned in a permuted form, according to the\n", + " optional return value P. For example, given the matrix 'a = [1, 2;\n", + " 3, 4]',\n", + "\n", + " [l, u, p] = lu (A)\n", + "\n", + " returns\n", + "\n", + " l =\n", + "\n", + " 1.00000 0.00000\n", + " 0.33333 1.00000\n", + "\n", + " u =\n", + "\n", + " 3.00000 4.00000\n", + " 0.00000 0.66667\n", + "\n", + " p =\n", + "\n", + " 0 1\n", + " 1 0\n", + "\n", + " The matrix is not required to be square.\n", + "\n", + " When called with two or three output arguments and a spare input\n", + " matrix, 'lu' does not attempt to perform sparsity preserving column\n", + " permutations. Called with a fourth output argument, the sparsity\n", + " preserving column transformation Q is returned, such that 'P * A *\n", + " Q = L * U'.\n", + "\n", + " Called with a fifth output argument and a sparse input matrix, 'lu'\n", + " attempts to use a scaling factor R on the input matrix such that 'P\n", + " * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n", + " stable factorization.\n", + "\n", + " An additional input argument THRES, that defines the pivoting\n", + " threshold can be given. THRES can be a scalar, in which case it\n", + " defines the UMFPACK pivoting tolerance for both symmetric and\n", + " unsymmetric cases. If THRES is a 2-element vector, then the first\n", + " element defines the pivoting tolerance for the unsymmetric UMFPACK\n", + " pivoting strategy and the second for the symmetric strategy. By\n", + " default, the values defined by 'spparms' are used ([0.1, 0.001]).\n", + "\n", + " Given the string argument \"vector\", 'lu' returns the values of P\n", + " and Q as vector values, such that for full matrix, 'A (P,:) = L *\n", + " U', and 'R(P,:) * A (:, Q) = L * U'.\n", + "\n", + " With two output arguments, returns the permuted forms of the upper\n", + " and lower triangular matrices, such that 'A = L * U'. With one\n", + " output argument Y, then the matrix returned by the LAPACK routines\n", + " is returned. If the input matrix is sparse then the matrix L is\n", + " embedded into U to give a return value similar to the full case.\n", + " For both full and sparse matrices, 'lu' loses the permutation\n", + " information.\n", + "\n", + " See also: luupdate, ilu, chol, hess, qr, qz, schur, svd.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help lu" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tLU decomp\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tLU decomp\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tOctave \\\\\n", + "\n", + "\t\n", + "\t\tOctave \\\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% time LU solution vs backslash\n", + "t_lu=zeros(100,1);\n", + "t_bs=zeros(100,1);\n", + "for N=1:100\n", + " A=rand(N,N);\n", + " y=rand(N,1);\n", + " [L,U,P]=lu(A);\n", + "\n", + " tic; d=inv(L)*y; x=inv(U)*d; t_lu(N)=toc;\n", + "\n", + " tic; x=inv(A)*y; t_bs(N)=toc;\n", + "end\n", + "plot([1:100],t_lu,[1:100],t_bs) \n", + "legend('LU decomp','Octave \\\\')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n", + "\n", + "![Springs-masses](../lecture_09/mass_springs.svg)\n", + "\n", + "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n", + "\n", + "$m_{1}g+k(x_{2}-x_{1})-kx_{1}=0$\n", + "\n", + "$m_{2}g+k(x_{3}-x_{2})-k(x_{2}-x_{1})=0$\n", + "\n", + "$m_{3}g+k(x_{4}-x_{3})-k(x_{3}-x_{2})=0$\n", + "\n", + "$m_{4}g-k(x_{4}-x_{3})=0$\n", + "\n", + "in matrix form:\n", + "\n", + "$\\left[ \\begin{array}{cccc}\n", + "2k & -k & 0 & 0 \\\\\n", + "-k & 2k & -k & 0 \\\\\n", + "0 & -k & 2k & -k \\\\\n", + "0 & 0 & -k & k \\end{array} \\right]\n", + "\\left[ \\begin{array}{c}\n", + "x_{1} \\\\\n", + "x_{2} \\\\\n", + "x_{3} \\\\\n", + "x_{4} \\end{array} \\right]=\n", + "\\left[ \\begin{array}{c}\n", + "m_{1}g \\\\\n", + "m_{2}g \\\\\n", + "m_{3}g \\\\\n", + "m_{4}g \\end{array} \\right]$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 20 -10 0 0\n", + " -10 20 -10 0\n", + " 0 -10 20 -10\n", + " 0 0 -10 10\n", + "\n", + "y =\n", + "\n", + " 9.8100\n", + " 19.6200\n", + " 29.4300\n", + " 39.2400\n", + "\n" + ] + } + ], + "source": [ + "k=10; % N/m\n", + "m1=1; % kg\n", + "m2=2;\n", + "m3=3;\n", + "m4=4;\n", + "g=9.81; % m/s^2\n", + "K=[2*k -k 0 0; -k 2*k -k 0; 0 -k 2*k -k; 0 0 -k k]\n", + "y=[m1*g;m2*g;m3*g;m4*g]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This matrix, K, is symmetric. \n", + "\n", + "`K(i,j)==K(j,i)`\n", + "\n", + "Now we can use,\n", + "\n", + "## Cholesky Factorization\n", + "\n", + "We can decompose the matrix, K into two matrices, $U$ and $U^{T}$, where\n", + "\n", + "$K=U^{T}U$\n", + "\n", + "each of the components of U can be calculated with the following equations:\n", + "\n", + "$u_{ii}=\\sqrt{a_{ii}-\\sum_{k=1}^{i-1}u_{ki}^{2}}$\n", + "\n", + "$u_{ij}=\\frac{a_{ij}-\\sum_{k=1}^{i-1}u_{ki}u_{kj}}{u_{ii}}$\n", + "\n", + "so for K" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 20 -10 0 0\n", + " -10 20 -10 0\n", + " 0 -10 20 -10\n", + " 0 0 -10 10\n", + "\n" + ] + } + ], + "source": [ + "K" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "u11 = 4.4721\n", + "u12 = -2.2361\n", + "u13 = 0\n", + "u14 = 0\n", + "u22 = 3.8730\n", + "u23 = -2.5820\n", + "u24 = 0\n", + "u33 = 3.6515\n", + "u34 = -2.7386\n", + "u44 = 1.5811\n", + "U =\n", + "\n", + " 4.47214 -2.23607 0.00000 0.00000\n", + " 0.00000 3.87298 -2.58199 0.00000\n", + " 0.00000 0.00000 3.65148 -2.73861\n", + " 0.00000 0.00000 0.00000 1.58114\n", + "\n" + ] + } + ], + "source": [ + "u11=sqrt(K(1,1))\n", + "u12=(K(1,2))/u11\n", + "u13=(K(1,3))/u11\n", + "u14=(K(1,4))/u11\n", + "u22=sqrt(K(2,2)-u12^2)\n", + "u23=(K(2,3)-u12*u13)/u22\n", + "u24=(K(2,4)-u12*u14)/u22\n", + "u33=sqrt(K(3,3)-u13^2-u23^2)\n", + "u34=(K(3,4)-u13*u14-u23*u24)/u33\n", + "u44=sqrt(K(4,4)-u14^2-u24^2-u34^2)\n", + "U=[u11,u12,u13,u14;0,u22,u23,u24;0,0,u33,u34;0,0,0,u44]" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans =\n", + "\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + "\n" + ] + } + ], + "source": [ + "(U'*U)'==U'*U" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "average time spent for Cholesky factored solution = 1.465964e-05+/-9.806001e-06\n", + "average time spent for backslash solution = 1.555967e-05+/-1.048561e-05\n" + ] + } + ], + "source": [ + "% time solution for Cholesky vs backslash\n", + "t_chol=zeros(1000,1);\n", + "t_bs=zeros(1000,1);\n", + "for i=1:1000\n", + " tic; d=U'*y; x=U\\d; t_chol(i)=toc;\n", + " tic; x=K\\y; t_bs(i)=toc;\n", + "end\n", + "fprintf('average time spent for Cholesky factored solution = %e+/-%e',mean(t_chol),std(t_chol))\n", + "\n", + "fprintf('average time spent for backslash solution = %e+/-%e',mean(t_bs),std(t_bs))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_12/LU_naive.m b/lecture_12/LU_naive.m new file mode 100644 index 0000000..92efde6 --- /dev/null +++ b/lecture_12/LU_naive.m @@ -0,0 +1,27 @@ +function [L, U] = LU_naive(A) +% GaussNaive: naive Gauss elimination +% x = GaussNaive(A,b): Gauss elimination without pivoting. +% input: +% A = coefficient matrix +% y = right hand side vector +% output: +% x = solution vector +[m,n] = size(A); +if m~=n, error('Matrix A must be square'); end +nb = n; +L=diag(ones(n,1)); +U=A; +% forward elimination +for k = 1:n-1 + for i = k+1:n + fik = U(i,k)/U(k,k); + L(i,k)=fik; + U(i,k:nb) = U(i,k:nb)-fik*U(k,k:nb); + end +end +%% back substitution +%x = zeros(n,1); +%x(n) = Aug(n,nb)/Aug(n,n); +%for i = n-1:-1:1 +% x(i) = (Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i); +%end diff --git a/lecture_12/chol_pre.png b/lecture_12/chol_pre.png new file mode 100644 index 0000000..67d2493 Binary files /dev/null and b/lecture_12/chol_pre.png differ diff --git a/lecture_12/det_L.png b/lecture_12/det_L.png new file mode 100644 index 0000000..aaac5a2 Binary files /dev/null and b/lecture_12/det_L.png differ diff --git a/lecture_12/lecture_12.aux b/lecture_12/lecture_12.aux new file mode 100644 index 0000000..c73dc14 --- /dev/null +++ b/lecture_12/lecture_12.aux @@ -0,0 +1,27 @@ +\relax +\providecommand\hyper@newdestlabel[2]{} +\providecommand\HyperFirstAtBeginDocument{\AtBeginDocument} +\HyperFirstAtBeginDocument{\ifx\hyper@anchor\@undefined +\global\let\oldcontentsline\contentsline +\gdef\contentsline#1#2#3#4{\oldcontentsline{#1}{#2}{#3}} +\global\let\oldnewlabel\newlabel +\gdef\newlabel#1#2{\newlabelxx{#1}#2} +\gdef\newlabelxx#1#2#3#4#5#6{\oldnewlabel{#1}{{#2}{#3}}} +\AtEndDocument{\ifx\hyper@anchor\@undefined +\let\contentsline\oldcontentsline +\let\newlabel\oldnewlabel +\fi} +\fi} +\global\let\hyper@last\relax +\gdef\HyperFirstAtBeginDocument#1{#1} +\providecommand\HyField@AuxAddToFields[1]{} +\providecommand\HyField@AuxAddToCoFields[2]{} +\providecommand \oddpage@label [2]{} +\@writefile{toc}{\contentsline {subsection}{\numberline {0.1}My question from last class}{1}{subsection.0.1}} +\newlabel{my-question-from-last-class}{{0.1}{1}{My question from last class}{subsection.0.1}{}} +\@writefile{lof}{\contentsline {figure}{\numberline {1}{\ignorespaces q1\relax }}{1}{figure.caption.1}} +\@writefile{toc}{\contentsline {subsection}{\numberline {0.2}Your questions from last class}{1}{subsection.0.2}} +\newlabel{your-questions-from-last-class}{{0.2}{1}{Your questions from last class}{subsection.0.2}{}} +\@writefile{lof}{\contentsline {figure}{\numberline {2}{\ignorespaces q2\relax }}{2}{figure.caption.2}} +\@writefile{toc}{\contentsline {section}{\numberline {1}Matrix Inverse and Condition}{2}{section.1}} +\newlabel{matrix-inverse-and-condition}{{1}{2}{Matrix Inverse and Condition}{section.1}{}} diff --git a/lecture_12/lecture_12.ipynb b/lecture_12/lecture_12.ipynb new file mode 100644 index 0000000..dc44a21 --- /dev/null +++ b/lecture_12/lecture_12.ipynb @@ -0,0 +1,1123 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## My question from last class \n", + "\n", + "![q1](det_L.png)\n", + "\n", + "![q2](chol_pre.png)\n", + "\n", + "\n", + "## Your questions from last class\n", + "\n", + "1. Will the exam be more theoretical or problem based?\n", + "\n", + "2. Writing code is difficult \n", + "\n", + "3. What format can we expect for the midterm? \n", + "\n", + "2. Could we go over some example questions for the exam?\n", + "\n", + "3. Will the use of GitHub be tested on the Midterm exam? Or is it more focused on linear algebra techniques/what was covered in the lectures?\n", + "\n", + "4. This is not my strong suit, getting a bit overwhelmed with matrix multiplication.\n", + "\n", + "5. I forgot how much I learned in linear algebra.\n", + "\n", + "6. What's the most exciting project you've ever worked on with Matlab/Octave?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Matrix Inverse and Condition\n", + "\n", + "\n", + "Considering the same solution set:\n", + "\n", + "$y=Ax$\n", + "\n", + "If we know that $A^{-1}A=I$, then \n", + "\n", + "$A^{-1}y=A^{-1}Ax=x$\n", + "\n", + "so \n", + "\n", + "$x=A^{-1}y$\n", + "\n", + "Where, $A^{-1}$ is the inverse of matrix $A$.\n", + "\n", + "$2x_{1}+x_{2}=1$\n", + "\n", + "$x_{1}+3x_{2}=1$\n", + "\n", + "$Ax=y$\n", + "\n", + "$\\left[ \\begin{array}{cc}\n", + "2 & 1 \\\\\n", + "1 & 3 \\end{array} \\right]\n", + "\\left[\\begin{array}{c} \n", + "x_{1} \\\\ \n", + "x_{2} \\end{array}\\right]=\n", + "\\left[\\begin{array}{c} \n", + "1 \\\\\n", + "1\\end{array}\\right]$\n", + "\n", + "$A^{-1}=\\frac{1}{2*3-1*1}\\left[ \\begin{array}{cc}\n", + "3 & 1 \\\\\n", + "-1 & 2 \\end{array} \\right]=\n", + "\\left[ \\begin{array}{cc}\n", + "3/5 & -1/5 \\\\\n", + "-1/5 & 2/5 \\end{array} \\right]$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "invA =\n", + "\n", + " 0.60000 -0.20000\n", + " -0.20000 0.40000\n", + "\n", + "ans =\n", + "\n", + " 1.00000 0.00000\n", + " 0.00000 1.00000\n", + "\n", + "ans =\n", + "\n", + " 1.00000 0.00000\n", + " 0.00000 1.00000\n", + "\n" + ] + } + ], + "source": [ + "A=[2,1;1,3]\n", + "invA=1/5*[3,-1;-1,2]\n", + "\n", + "A*invA\n", + "invA*A" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "How did we know the inverse of A? \n", + "\n", + "for 2$\\times$2 matrices, it is always:\n", + "\n", + "$A=\\left[ \\begin{array}{cc}\n", + "A_{11} & A_{12} \\\\\n", + "A_{21} & A_{22} \\end{array} \\right]$\n", + "\n", + "$A^{-1}=\\frac{1}{det(A)}\\left[ \\begin{array}{cc}\n", + "A_{22} & -A_{12} \\\\\n", + "-A_{21} & A_{11} \\end{array} \\right]$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$AA^{-1}=\\frac{1}{A_{11}A_{22}-A_{21}A_{12}}\\left[ \\begin{array}{cc}\n", + "A_{11}A_{22}-A_{21}A_{12} & -A_{11}A_{12}+A_{12}A_{11} \\\\\n", + "A_{21}A_{22}-A_{22}A_{21} & -A_{21}A_{12}+A_{22}A_{11} \\end{array} \\right]\n", + "=\\left[ \\begin{array}{cc}\n", + "1 & 0 \\\\\n", + "0 & 1 \\end{array} \\right]$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What about bigger matrices?\n", + "\n", + "We can use the LU-decomposition\n", + "\n", + "$A=LU$\n", + "\n", + "$A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$\n", + "\n", + "if we divide $A^{-1}$ into n-column vectors, $a_{n}$, then\n", + "\n", + "$Aa_{1}=\\left[\\begin{array}{c} \n", + "1 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$Aa_{2}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "1 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$Aa_{n}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "1 \\end{array} \\right]$\n", + "\n", + "\n", + "Which we can solve for each $a_{n}$ with LU-decomposition, knowing the lower and upper triangular decompositions, then \n", + "\n", + "$A^{-1}=\\left[ \\begin{array}{cccc}\n", + "| & | & & | \\\\\n", + "a_{1} & a_{2} & \\cdots & a_{3} \\\\\n", + "| & | & & | \\end{array} \\right]$\n", + "\n", + "\n", + "$Ld_{1}=\\left[\\begin{array}{c} \n", + "1 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$;~Ua_{1}=d_{1}$\n", + "\n", + "$Ld_{2}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "1 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$;~Ua_{2}=d_{2}$\n", + "\n", + "$Ld_{n}=\\left[\\begin{array}{c} \n", + "0 \\\\ \n", + "1 \\\\ \n", + "\\vdots \\\\\n", + "n \\end{array} \\right]$\n", + "$;~Ua_{n}=d_{n}$\n", + "\n", + "Consider the following matrix:\n", + "\n", + "$A=\\left[ \\begin{array}{ccc}\n", + "2 & -1 & 0\\\\\n", + "-1 & 2 & -1\\\\\n", + "0 & -1 & 1 \\end{array} \\right]$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 -1 0\n", + " -1 2 -1\n", + " 0 -1 1\n", + "\n", + "U =\n", + "\n", + " 2.00000 -1.00000 0.00000\n", + " 0.00000 1.50000 -1.00000\n", + " 0.00000 -1.00000 1.00000\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " -0.50000 1.00000 0.00000\n", + " 0.00000 0.00000 1.00000\n", + "\n" + ] + } + ], + "source": [ + "A=[2,-1,0;-1,2,-1;0,-1,1]\n", + "U=A;\n", + "L=eye(3,3);\n", + "U(2,:)=U(2,:)-U(2,1)/U(1,1)*U(1,:)\n", + "L(2,1)=A(2,1)/A(1,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "L =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " -0.50000 1.00000 0.00000\n", + " 0.00000 -0.66667 1.00000\n", + "\n", + "U =\n", + "\n", + " 2.00000 -1.00000 0.00000\n", + " 0.00000 1.50000 -1.00000\n", + " 0.00000 0.00000 0.33333\n", + "\n" + ] + } + ], + "source": [ + "L(3,2)=U(3,2)/U(2,2)\n", + "U(3,:)=U(3,:)-U(3,2)/U(2,2)*U(2,:)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now solve for $d_1$ then $a_1$, $d_2$ then $a_2$, and $d_3$ then $a_{3}$\n", + "\n", + "$Ld_{1}=\\left[\\begin{array}{c} \n", + "1 \\\\ \n", + "0 \\\\ \n", + "\\vdots \\\\\n", + "0 \\end{array} \\right]$\n", + "$;~Ua_{1}=d_{1}$" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "d1 =\n", + "\n", + " 1.00000\n", + " 0.50000\n", + " 0.33333\n", + "\n" + ] + } + ], + "source": [ + "d1=zeros(3,1);\n", + "d1(1)=1;\n", + "d1(2)=0-L(2,1)*d1(1);\n", + "d1(3)=0-L(3,1)*d1(1)-L(3,2)*d1(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a1 =\n", + "\n", + " 1.00000\n", + " 1.00000\n", + " 1.00000\n", + "\n" + ] + } + ], + "source": [ + "a1=zeros(3,1);\n", + "a1(3)=d1(3)/U(3,3);\n", + "a1(2)=1/U(2,2)*(d1(2)-U(2,3)*a1(3));\n", + "a1(1)=1/U(1,1)*(d1(1)-U(1,2)*a1(2)-U(1,3)*a1(3))" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "d2 =\n", + "\n", + " 0.00000\n", + " 1.00000\n", + " 0.66667\n", + "\n" + ] + } + ], + "source": [ + "d2=zeros(3,1);\n", + "d2(1)=0;\n", + "d2(2)=1-L(2,1)*d2(1);\n", + "d2(3)=0-L(3,1)*d2(1)-L(3,2)*d2(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a2 =\n", + "\n", + " 1.0000\n", + " 2.0000\n", + " 2.0000\n", + "\n" + ] + } + ], + "source": [ + "a2=zeros(3,1);\n", + "a2(3)=d2(3)/U(3,3);\n", + "a2(2)=1/U(2,2)*(d2(2)-U(2,3)*a2(3));\n", + "a2(1)=1/U(1,1)*(d2(1)-U(1,2)*a2(2)-U(1,3)*a2(3))" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "d3 =\n", + "\n", + " 0\n", + " 0\n", + " 1\n", + "\n" + ] + } + ], + "source": [ + "d3=zeros(3,1);\n", + "d3(1)=0;\n", + "d3(2)=0-L(2,1)*d3(1);\n", + "d3(3)=1-L(3,1)*d3(1)-L(3,2)*d3(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "a3 =\n", + "\n", + " 1.00000\n", + " 2.00000\n", + " 3.00000\n", + "\n" + ] + } + ], + "source": [ + "a3=zeros(3,1);\n", + "a3(3)=d3(3)/U(3,3);\n", + "a3(2)=1/U(2,2)*(d3(2)-U(2,3)*a3(3));\n", + "a3(1)=1/U(1,1)*(d3(1)-U(1,2)*a3(2)-U(1,3)*a3(3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Final solution for $A^{-1}$ is $[a_{1}~a_{2}~a_{3}]$" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "invA =\n", + "\n", + " 1.00000 1.00000 1.00000\n", + " 1.00000 2.00000 2.00000\n", + " 1.00000 2.00000 3.00000\n", + "\n", + "ans =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " 0.00000 1.00000 -0.00000\n", + " -0.00000 -0.00000 1.00000\n", + "\n" + ] + } + ], + "source": [ + "invA=[a1,a2,a3]\n", + "A*invA" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now the solution of $x$ to $Ax=y$ is $x=A^{-1}y$" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "y =\n", + "\n", + " 1\n", + " 2\n", + " 3\n", + "\n", + "x =\n", + "\n", + " 6.0000\n", + " 11.0000\n", + " 14.0000\n", + "\n", + "xbs =\n", + "\n", + " 6.0000\n", + " 11.0000\n", + " 14.0000\n", + "\n", + "ans =\n", + "\n", + " -3.5527e-15\n", + " -8.8818e-15\n", + " -1.0658e-14\n", + "\n", + "ans = 2.2204e-16\n" + ] + } + ], + "source": [ + "y=[1;2;3]\n", + "x=invA*y\n", + "xbs=A\\y\n", + "x-xbs\n", + "eps" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tinversion\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tinversion\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tbackslash\n", + "\n", + "\t\n", + "\t\tbackslash\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tmultiplication\n", + "\n", + "\t\n", + "\t\tmultiplication\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=100;\n", + "n=[1:N];\n", + "t_inv=zeros(N,1);\n", + "t_bs=zeros(N,1);\n", + "t_mult=zeros(N,1);\n", + "for i=1:N\n", + " A=rand(i,i);\n", + " tic\n", + " invA=inv(A);\n", + " t_inv(i)=toc;\n", + " b=rand(i,1);\n", + " tic;\n", + " x=A\\b;\n", + " t_bs(i)=toc;\n", + " tic;\n", + " x=invA*b;\n", + " t_mult(i)=toc;\n", + "end\n", + "plot(n,t_inv,n,t_bs,n,t_mult)\n", + "axis([0 100 0 0.002])\n", + "legend('inversion','backslash','multiplication','Location','NorthWest')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Condition of a matrix \n", + "### *just checked in to see what condition my condition was in*\n", + "### Matrix norms\n", + "\n", + "The Euclidean norm of a vector is measure of the magnitude (in 3D this would be: $|x|=\\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$) in general the equation is:\n", + "\n", + "$||x||_{e}=\\sqrt{\\sum_{i=1}^{n}x_{i}^{2}}$\n", + "\n", + "For a matrix, A, the same norm is called the Frobenius norm:\n", + "\n", + "$||A||_{f}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{2}}$\n", + "\n", + "In general we can calculate any $p$-norm where\n", + "\n", + "$||A||_{p}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{p}}$\n", + "\n", + "so the p=1, 1-norm is \n", + "\n", + "$||A||_{1}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{1}}=\\sum_{i=1}^{n}\\sum_{i=1}^{m}|A_{i,j}|$\n", + "\n", + "$||A||_{\\infty}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{\\infty}}=\\max_{1\\le i \\le n}\\sum_{j=1}^{m}|A_{i,j}|$\n", + "\n", + "### Condition of Matrix\n", + "\n", + "The matrix condition is the product of \n", + "\n", + "$Cond(A) = ||A||\\cdot||A^{-1}||$ \n", + "\n", + "So each norm will have a different condition number, but the limit is $Cond(A)\\ge 1$\n", + "\n", + "An estimate of the rounding error is based on the condition of A:\n", + "\n", + "$\\frac{||\\Delta x||}{x} \\le Cond(A) \\frac{||\\Delta A||}{||A||}$\n", + "\n", + "So if the coefficients of A have accuracy to $10^{-t}\n", + "\n", + "and the condition of A, $Cond(A)=10^{c}$\n", + "\n", + "then the solution for x can have rounding errors up to $10^{c-t}$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 1.00000 0.50000 0.33333\n", + " 0.50000 0.33333 0.25000\n", + " 0.33333 0.25000 0.20000\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000 0.00000\n", + " 0.50000 1.00000 0.00000\n", + " 0.33333 1.00000 1.00000\n", + "\n", + "U =\n", + "\n", + " 1.00000 0.50000 0.33333\n", + " 0.00000 0.08333 0.08333\n", + " 0.00000 -0.00000 0.00556\n", + "\n" + ] + } + ], + "source": [ + "A=[1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5]\n", + "[L,U]=LU_naive(A)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "Then, $A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$\n", + "\n", + "$Ld_{1}=\\left[\\begin{array}{c}\n", + "1 \\\\\n", + "0 \\\\\n", + "0 \\end{array}\\right]$, $Ux_{1}=d_{1}$ ..." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "invA =\n", + "\n", + " 9.0000 -36.0000 30.0000\n", + " -36.0000 192.0000 -180.0000\n", + " 30.0000 -180.0000 180.0000\n", + "\n", + "ans =\n", + "\n", + " 1.0000e+00 3.5527e-15 2.9976e-15\n", + " -1.3249e-14 1.0000e+00 -9.1038e-15\n", + " 8.5117e-15 7.1054e-15 1.0000e+00\n", + "\n" + ] + } + ], + "source": [ + "invA=zeros(3,3);\n", + "d1=L\\[1;0;0];\n", + "d2=L\\[0;1;0];\n", + "d3=L\\[0;0;1];\n", + "invA(:,1)=U\\d1;\n", + "invA(:,2)=U\\d2;\n", + "invA(:,3)=U\\d3\n", + "invA*A" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Find the condition of A, $cond(A)$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "normf_A = 1.4136\n", + "normf_invA = 372.21\n", + "cond_f_A = 526.16\n", + "ans = 1.4136\n", + "norm1_A = 1.8333\n", + "norm1_invA = 30.000\n", + "ans = 1.8333\n", + "cond_1_A = 55.000\n", + "norminf_A = 1.8333\n", + "norminf_invA = 30.000\n", + "ans = 1.8333\n", + "cond_inf_A = 55.000\n" + ] + } + ], + "source": [ + "% Frobenius norm\n", + "normf_A = sqrt(sum(sum(A.^2)))\n", + "normf_invA = sqrt(sum(sum(invA.^2)))\n", + "\n", + "cond_f_A = normf_A*normf_invA\n", + "\n", + "norm(A,'fro')\n", + "\n", + "% p=1, column sum norm\n", + "norm1_A = max(sum(A,2))\n", + "norm1_invA = max(sum(invA,2))\n", + "norm(A,1)\n", + "\n", + "cond_1_A=norm1_A*norm1_invA\n", + "\n", + "% p=inf, row sum norm\n", + "norminf_A = max(sum(A,1))\n", + "norminf_invA = max(sum(invA,1))\n", + "norm(A,inf)\n", + "\n", + "cond_inf_A=norminf_A*norminf_invA\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with spring constants $K_{i}$. What does a high condition number mean for this problem? \n", + "\n", + "![Springs-masses](../lecture_09/mass_springs.png)\n", + "\n", + "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n", + "\n", + "$m_{1}g+k_{2}(x_{2}-x_{1})-k_{1}x_{1}=0$\n", + "\n", + "$m_{2}g+k_{3}(x_{3}-x_{2})-k_{2}(x_{2}-x_{1})=0$\n", + "\n", + "$m_{3}g+k_{4}(x_{4}-x_{3})-k_{3}(x_{3}-x_{2})=0$\n", + "\n", + "$m_{4}g-k_{4}(x_{4}-x_{3})=0$\n", + "\n", + "in matrix form:\n", + "\n", + "$\\left[ \\begin{array}{cccc}\n", + "k_{1}+k_{2} & -k_{2} & 0 & 0 \\\\\n", + "-k_{2} & k_{2}+k_{3} & -k_{3} & 0 \\\\\n", + "0 & -k_{3} & k_{3}+k_{4} & -k_{4} \\\\\n", + "0 & 0 & -k_{4} & k_{4} \\end{array} \\right]\n", + "\\left[ \\begin{array}{c}\n", + "x_{1} \\\\\n", + "x_{2} \\\\\n", + "x_{3} \\\\\n", + "x_{4} \\end{array} \\right]=\n", + "\\left[ \\begin{array}{c}\n", + "m_{1}g \\\\\n", + "m_{2}g \\\\\n", + "m_{3}g \\\\\n", + "m_{4}g \\end{array} \\right]$" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 100010 -100000 0 0\n", + " -100000 100010 -10 0\n", + " 0 -10 11 -1\n", + " 0 0 -1 1\n", + "\n", + "y =\n", + "\n", + " 9.8100\n", + " 19.6200\n", + " 29.4300\n", + " 39.2400\n", + "\n" + ] + } + ], + "source": [ + "k1=10; % N/m\n", + "k2=100000;\n", + "k3=10;\n", + "k4=1;\n", + "m1=1; % kg\n", + "m2=2;\n", + "m3=3;\n", + "m4=4;\n", + "g=9.81; % m/s^2\n", + "K=[k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4]\n", + "y=[m1*g;m2*g;m3*g;m4*g]" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 3.2004e+05\n", + "ans = 3.2004e+05\n", + "ans = 2.5925e+05\n", + "ans = 2.5293e+05\n" + ] + } + ], + "source": [ + "cond(K,inf)\n", + "cond(K,1)\n", + "cond(K,'fro')\n", + "cond(K,2)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "e =\n", + "\n", + " 7.9078e-01\n", + " 3.5881e+00\n", + " 1.7621e+01\n", + " 2.0001e+05\n", + "\n", + "ans = 2.5293e+05\n" + ] + } + ], + "source": [ + "e=eig(K)\n", + "max(e)/min(e)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_12/lecture_12.md b/lecture_12/lecture_12.md new file mode 100644 index 0000000..9befcbc --- /dev/null +++ b/lecture_12/lecture_12.md @@ -0,0 +1,674 @@ + + +```octave +%plot --format svg +``` + + +```octave +setdefaults +``` + +## My question from last class + +![q1](det_L.png) + +![q2](chol_pre.png) + + +## Your questions from last class + +1. Will the exam be more theoretical or problem based? + +2. Writing code is difficult + +3. What format can we expect for the midterm? + +2. Could we go over some example questions for the exam? + +3. Will the use of GitHub be tested on the Midterm exam? Or is it more focused on linear algebra techniques/what was covered in the lectures? + +4. This is not my strong suit, getting a bit overwhelmed with matrix multiplication. + +5. I forgot how much I learned in linear algebra. + +6. What's the most exciting project you've ever worked on with Matlab/Octave? + +# Matrix Inverse and Condition + + +Considering the same solution set: + +$y=Ax$ + +If we know that $A^{-1}A=I$, then + +$A^{-1}y=A^{-1}Ax=x$ + +so + +$x=A^{-1}y$ + +Where, $A^{-1}$ is the inverse of matrix $A$. + +$2x_{1}+x_{2}=1$ + +$x_{1}+3x_{2}=1$ + +$Ax=y$ + +$\left[ \begin{array}{cc} +2 & 1 \\ +1 & 3 \end{array} \right] +\left[\begin{array}{c} +x_{1} \\ +x_{2} \end{array}\right]= +\left[\begin{array}{c} +1 \\ +1\end{array}\right]$ + +$A^{-1}=\frac{1}{2*3-1*1}\left[ \begin{array}{cc} +3 & 1 \\ +-1 & 2 \end{array} \right]= +\left[ \begin{array}{cc} +3/5 & -1/5 \\ +-1/5 & 2/5 \end{array} \right]$ + + + +```octave +A=[2,1;1,3] +invA=1/5*[3,-1;-1,2] + +A*invA +invA*A +``` + + A = + + 2 1 + 1 3 + + invA = + + 0.60000 -0.20000 + -0.20000 0.40000 + + ans = + + 1.00000 0.00000 + 0.00000 1.00000 + + ans = + + 1.00000 0.00000 + 0.00000 1.00000 + + + +How did we know the inverse of A? + +for 2$\times$2 matrices, it is always: + +$A=\left[ \begin{array}{cc} +A_{11} & A_{12} \\ +A_{21} & A_{22} \end{array} \right]$ + +$A^{-1}=\frac{1}{det(A)}\left[ \begin{array}{cc} +A_{22} & -A_{12} \\ +-A_{21} & A_{11} \end{array} \right]$ + +$AA^{-1}=\frac{1}{A_{11}A_{22}-A_{21}A_{12}}\left[ \begin{array}{cc} +A_{11}A_{22}-A_{21}A_{12} & -A_{11}A_{12}+A_{12}A_{11} \\ +A_{21}A_{22}-A_{22}A_{21} & -A_{21}A_{12}+A_{22}A_{11} \end{array} \right] +=\left[ \begin{array}{cc} +1 & 0 \\ +0 & 1 \end{array} \right]$ + +What about bigger matrices? + +We can use the LU-decomposition + +$A=LU$ + +$A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$ + +if we divide $A^{-1}$ into n-column vectors, $a_{n}$, then + +$Aa_{1}=\left[\begin{array}{c} +1 \\ +0 \\ +\vdots \\ +0 \end{array} \right]$ +$Aa_{2}=\left[\begin{array}{c} +0 \\ +1 \\ +\vdots \\ +0 \end{array} \right]$ +$Aa_{n}=\left[\begin{array}{c} +0 \\ +0 \\ +\vdots \\ +1 \end{array} \right]$ + + +Which we can solve for each $a_{n}$ with LU-decomposition, knowing the lower and upper triangular decompositions, then + +$A^{-1}=\left[ \begin{array}{cccc} +| & | & & | \\ +a_{1} & a_{2} & \cdots & a_{3} \\ +| & | & & | \end{array} \right]$ + + +$Ld_{1}=\left[\begin{array}{c} +1 \\ +0 \\ +\vdots \\ +0 \end{array} \right]$ +$;~Ua_{1}=d_{1}$ + +$Ld_{2}=\left[\begin{array}{c} +0 \\ +1 \\ +\vdots \\ +0 \end{array} \right]$ +$;~Ua_{2}=d_{2}$ + +$Ld_{n}=\left[\begin{array}{c} +0 \\ +1 \\ +\vdots \\ +n \end{array} \right]$ +$;~Ua_{n}=d_{n}$ + +Consider the following matrix: + +$A=\left[ \begin{array}{ccc} +2 & -1 & 0\\ +-1 & 2 & -1\\ +0 & -1 & 1 \end{array} \right]$ + + + +```octave +A=[2,-1,0;-1,2,-1;0,-1,1] +U=A; +L=eye(3,3); +U(2,:)=U(2,:)-U(2,1)/U(1,1)*U(1,:) +L(2,1)=A(2,1)/A(1,1) +``` + + A = + + 2 -1 0 + -1 2 -1 + 0 -1 1 + + U = + + 2.00000 -1.00000 0.00000 + 0.00000 1.50000 -1.00000 + 0.00000 -1.00000 1.00000 + + L = + + 1.00000 0.00000 0.00000 + -0.50000 1.00000 0.00000 + 0.00000 0.00000 1.00000 + + + + +```octave +L(3,2)=U(3,2)/U(2,2) +U(3,:)=U(3,:)-U(3,2)/U(2,2)*U(2,:) + +``` + + L = + + 1.00000 0.00000 0.00000 + -0.50000 1.00000 0.00000 + 0.00000 -0.66667 1.00000 + + U = + + 2.00000 -1.00000 0.00000 + 0.00000 1.50000 -1.00000 + 0.00000 0.00000 0.33333 + + + +Now solve for $d_1$ then $a_1$, $d_2$ then $a_2$, and $d_3$ then $a_{3}$ + +$Ld_{1}=\left[\begin{array}{c} +1 \\ +0 \\ +\vdots \\ +0 \end{array} \right]$ +$;~Ua_{1}=d_{1}$ + + +```octave +d1=zeros(3,1); +d1(1)=1; +d1(2)=0-L(2,1)*d1(1); +d1(3)=0-L(3,1)*d1(1)-L(3,2)*d1(2) +``` + + d1 = + + 1.00000 + 0.50000 + 0.33333 + + + + +```octave +a1=zeros(3,1); +a1(3)=d1(3)/U(3,3); +a1(2)=1/U(2,2)*(d1(2)-U(2,3)*a1(3)); +a1(1)=1/U(1,1)*(d1(1)-U(1,2)*a1(2)-U(1,3)*a1(3)) +``` + + a1 = + + 1.00000 + 1.00000 + 1.00000 + + + + +```octave +d2=zeros(3,1); +d2(1)=0; +d2(2)=1-L(2,1)*d2(1); +d2(3)=0-L(3,1)*d2(1)-L(3,2)*d2(2) +``` + + d2 = + + 0.00000 + 1.00000 + 0.66667 + + + + +```octave +a2=zeros(3,1); +a2(3)=d2(3)/U(3,3); +a2(2)=1/U(2,2)*(d2(2)-U(2,3)*a2(3)); +a2(1)=1/U(1,1)*(d2(1)-U(1,2)*a2(2)-U(1,3)*a2(3)) +``` + + a2 = + + 1.0000 + 2.0000 + 2.0000 + + + + +```octave +d3=zeros(3,1); +d3(1)=0; +d3(2)=0-L(2,1)*d3(1); +d3(3)=1-L(3,1)*d3(1)-L(3,2)*d3(2) +``` + + d3 = + + 0 + 0 + 1 + + + + +```octave +a3=zeros(3,1); +a3(3)=d3(3)/U(3,3); +a3(2)=1/U(2,2)*(d3(2)-U(2,3)*a3(3)); +a3(1)=1/U(1,1)*(d3(1)-U(1,2)*a3(2)-U(1,3)*a3(3)) +``` + + a3 = + + 1.00000 + 2.00000 + 3.00000 + + + +Final solution for $A^{-1}$ is $[a_{1}~a_{2}~a_{3}]$ + + +```octave +invA=[a1,a2,a3] +A*invA +``` + + invA = + + 1.00000 1.00000 1.00000 + 1.00000 2.00000 2.00000 + 1.00000 2.00000 3.00000 + + ans = + + 1.00000 0.00000 0.00000 + 0.00000 1.00000 -0.00000 + -0.00000 -0.00000 1.00000 + + + +Now the solution of $x$ to $Ax=y$ is $x=A^{-1}y$ + + +```octave +y=[1;2;3] +x=invA*y +xbs=A\y +x-xbs +eps +``` + + y = + + 1 + 2 + 3 + + x = + + 6.0000 + 11.0000 + 14.0000 + + xbs = + + 6.0000 + 11.0000 + 14.0000 + + ans = + + -3.5527e-15 + -8.8818e-15 + -1.0658e-14 + + ans = 2.2204e-16 + + + +```octave +N=100; +n=[1:N]; +t_inv=zeros(N,1); +t_bs=zeros(N,1); +t_mult=zeros(N,1); +for i=1:N + A=rand(i,i); + tic + invA=inv(A); + t_inv(i)=toc; + b=rand(i,1); + tic; + x=A\b; + t_bs(i)=toc; + tic; + x=invA*b; + t_mult(i)=toc; +end +plot(n,t_inv,n,t_bs,n,t_mult) +axis([0 100 0 0.002]) +legend('inversion','backslash','multiplication','Location','NorthWest') +``` + + +![svg](lecture_12_files/lecture_12_21_0.svg) + + +## Condition of a matrix +### *just checked in to see what condition my condition was in* +### Matrix norms + +The Euclidean norm of a vector is measure of the magnitude (in 3D this would be: $|x|=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$) in general the equation is: + +$||x||_{e}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$ + +For a matrix, A, the same norm is called the Frobenius norm: + +$||A||_{f}=\sqrt{\sum_{i=1}^{n}\sum_{i=1}^{m}A_{i,j}^{2}}$ + +In general we can calculate any $p$-norm where + +$||A||_{p}=\sqrt{\sum_{i=1}^{n}\sum_{i=1}^{m}A_{i,j}^{p}}$ + +so the p=1, 1-norm is + +$||A||_{1}=\sqrt{\sum_{i=1}^{n}\sum_{i=1}^{m}A_{i,j}^{1}}=\sum_{i=1}^{n}\sum_{i=1}^{m}|A_{i,j}|$ + +$||A||_{\infty}=\sqrt{\sum_{i=1}^{n}\sum_{i=1}^{m}A_{i,j}^{\infty}}=\max_{1\le i \le n}\sum_{j=1}^{m}|A_{i,j}|$ + +### Condition of Matrix + +The matrix condition is the product of + +$Cond(A) = ||A||\cdot||A^{-1}||$ + +So each norm will have a different condition number, but the limit is $Cond(A)\ge 1$ + +An estimate of the rounding error is based on the condition of A: + +$\frac{||\Delta x||}{x} \le Cond(A) \frac{||\Delta A||}{||A||}$ + +So if the coefficients of A have accuracy to $10^{-t} + +and the condition of A, $Cond(A)=10^{c}$ + +then the solution for x can have rounding errors up to $10^{c-t}$ + + + +```octave +A=[1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5] +[L,U]=LU_naive(A) +``` + + A = + + 1.00000 0.50000 0.33333 + 0.50000 0.33333 0.25000 + 0.33333 0.25000 0.20000 + + L = + + 1.00000 0.00000 0.00000 + 0.50000 1.00000 0.00000 + 0.33333 1.00000 1.00000 + + U = + + 1.00000 0.50000 0.33333 + 0.00000 0.08333 0.08333 + 0.00000 -0.00000 0.00556 + + + +Then, $A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$ + +$Ld_{1}=\left[\begin{array}{c} +1 \\ +0 \\ +0 \end{array}\right]$, $Ux_{1}=d_{1}$ ... + + +```octave +invA=zeros(3,3); +d1=L\[1;0;0]; +d2=L\[0;1;0]; +d3=L\[0;0;1]; +invA(:,1)=U\d1; +invA(:,2)=U\d2; +invA(:,3)=U\d3 +invA*A +``` + + invA = + + 9.0000 -36.0000 30.0000 + -36.0000 192.0000 -180.0000 + 30.0000 -180.0000 180.0000 + + ans = + + 1.0000e+00 3.5527e-15 2.9976e-15 + -1.3249e-14 1.0000e+00 -9.1038e-15 + 8.5117e-15 7.1054e-15 1.0000e+00 + + + +Find the condition of A, $cond(A)$ + + +```octave +% Frobenius norm +normf_A = sqrt(sum(sum(A.^2))) +normf_invA = sqrt(sum(sum(invA.^2))) + +cond_f_A = normf_A*normf_invA + +norm(A,'fro') + +% p=1, column sum norm +norm1_A = max(sum(A,2)) +norm1_invA = max(sum(invA,2)) +norm(A,1) + +cond_1_A=norm1_A*norm1_invA + +% p=inf, row sum norm +norminf_A = max(sum(A,1)) +norminf_invA = max(sum(invA,1)) +norm(A,inf) + +cond_inf_A=norminf_A*norminf_invA + +``` + + normf_A = 1.4136 + normf_invA = 372.21 + cond_f_A = 526.16 + ans = 1.4136 + norm1_A = 1.8333 + norm1_invA = 30.000 + ans = 1.8333 + cond_1_A = 55.000 + norminf_A = 1.8333 + norminf_invA = 30.000 + ans = 1.8333 + cond_inf_A = 55.000 + + +Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with spring constants $K_{i}$. What does a high condition number mean for this problem? + +![Springs-masses](../lecture_09/mass_springs.png) + +The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass: + +$m_{1}g+k_{2}(x_{2}-x_{1})-k_{1}x_{1}=0$ + +$m_{2}g+k_{3}(x_{3}-x_{2})-k_{2}(x_{2}-x_{1})=0$ + +$m_{3}g+k_{4}(x_{4}-x_{3})-k_{3}(x_{3}-x_{2})=0$ + +$m_{4}g-k_{4}(x_{4}-x_{3})=0$ + +in matrix form: + +$\left[ \begin{array}{cccc} +k_{1}+k_{2} & -k_{2} & 0 & 0 \\ +-k_{2} & k_{2}+k_{3} & -k_{3} & 0 \\ +0 & -k_{3} & k_{3}+k_{4} & -k_{4} \\ +0 & 0 & -k_{4} & k_{4} \end{array} \right] +\left[ \begin{array}{c} +x_{1} \\ +x_{2} \\ +x_{3} \\ +x_{4} \end{array} \right]= +\left[ \begin{array}{c} +m_{1}g \\ +m_{2}g \\ +m_{3}g \\ +m_{4}g \end{array} \right]$ + + +```octave +k1=10; % N/m +k2=100000; +k3=10; +k4=1; +m1=1; % kg +m2=2; +m3=3; +m4=4; +g=9.81; % m/s^2 +K=[k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4] +y=[m1*g;m2*g;m3*g;m4*g] +``` + + K = + + 100010 -100000 0 0 + -100000 100010 -10 0 + 0 -10 11 -1 + 0 0 -1 1 + + y = + + 9.8100 + 19.6200 + 29.4300 + 39.2400 + + + + +```octave +cond(K,inf) +cond(K,1) +cond(K,'fro') +cond(K,2) +``` + + ans = 3.2004e+05 + ans = 3.2004e+05 + ans = 2.5925e+05 + ans = 2.5293e+05 + + + +```octave +e=eig(K) +max(e)/min(e) +``` + + e = + + 7.9078e-01 + 3.5881e+00 + 1.7621e+01 + 2.0001e+05 + + ans = 2.5293e+05 + + + +```octave + +``` diff --git a/lecture_12/lecture_12.pdf b/lecture_12/lecture_12.pdf new file mode 100644 index 0000000..ed0f569 Binary files /dev/null and b/lecture_12/lecture_12.pdf differ diff --git a/lecture_12/lecture_12_files/lecture_12_21_0.pdf b/lecture_12/lecture_12_files/lecture_12_21_0.pdf new file mode 100644 index 0000000..6480131 Binary files /dev/null and b/lecture_12/lecture_12_files/lecture_12_21_0.pdf differ diff --git a/lecture_12/lecture_12_files/lecture_12_21_0.svg b/lecture_12/lecture_12_files/lecture_12_21_0.svg new file mode 100644 index 0000000..44ef608 --- /dev/null +++ b/lecture_12/lecture_12_files/lecture_12_21_0.svg @@ -0,0 +1,148 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0 + + + + + 20 + + + + + 40 + + + + + 60 + + + + + 80 + + + + + 100 + + + + + + + + + + + + + inversion + + + + + inversion + + + + + + backslash + + + backslash + + + + + + multiplication + + + multiplication + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_12/nohup.out b/lecture_12/nohup.out new file mode 100644 index 0000000..f0d30e1 --- /dev/null +++ b/lecture_12/nohup.out @@ -0,0 +1,2 @@ + +(evince:8926): Gtk-WARNING **: gtk_widget_size_allocate(): attempt to allocate widget with width -71 and height 20 diff --git a/lecture_12/octave-workspace b/lecture_12/octave-workspace new file mode 100644 index 0000000..41ef164 Binary files /dev/null and b/lecture_12/octave-workspace differ diff --git a/lecture_13/.ipynb_checkpoints/lecture_12-checkpoint.ipynb b/lecture_13/.ipynb_checkpoints/lecture_12-checkpoint.ipynb new file mode 100644 index 0000000..6a5dc68 --- /dev/null +++ b/lecture_13/.ipynb_checkpoints/lecture_12-checkpoint.ipynb @@ -0,0 +1,737 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## My question from last class \n", + "\n", + "\n", + "## Your questions from last class\n", + "\n", + "1. Need more linear algebra review\n", + " \n", + " -We will keep doing Linear Algebra, try the practice problems in 'linear_algebra'\n", + "\n", + "2. How do I do HW3? \n", + " \n", + " -demo today\n", + "\n", + "3. For hw4 is the spring constant (K) suppose to be given? \n", + " \n", + " -yes, its 30 N/m\n", + " \n", + "4. Deapool or Joker?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# LU Decomposition\n", + "### efficient storage of matrices for solutions\n", + "\n", + "Considering the same solution set:\n", + "\n", + "$y=Ax$\n", + "\n", + "Assume that we can perform Gauss elimination and achieve this formula:\n", + "\n", + "$Ux=d$ \n", + "\n", + "Where, $U$ is an upper triangular matrix that we derived from Gauss elimination and $d$ is the set of dependent variables after Gauss elimination. \n", + "\n", + "Assume there is a lower triangular matrix, $L$, with ones on the diagonal and same dimensions of $U$ and the following is true:\n", + "\n", + "$L(Ux-d)=Ax-y=0$\n", + "\n", + "Now, $Ax=LUx$, so $A=LU$, and $y=Ld$.\n", + "\n", + "$2x_{1}+x_{2}=1$\n", + "\n", + "$x_{1}+3x_{2}=1$\n", + "\n", + "\n", + "$\\left[ \\begin{array}{cc}\n", + "2 & 1 \\\\\n", + "1 & 3 \\end{array} \\right]\n", + "\\left[\\begin{array}{c} \n", + "x_{1} \\\\ \n", + "x_{2} \\end{array}\\right]=\n", + "\\left[\\begin{array}{c} \n", + "1 \\\\\n", + "1\\end{array}\\right]$\n", + "\n", + "f21=0.5\n", + "\n", + "A(2,1)=1-1 = 0 \n", + "\n", + "A(2,2)=3-0.5=2.5\n", + "\n", + "y(2)=1-0.5=0.5\n", + "\n", + "$L(Ux-d)=\n", + "\\left[ \\begin{array}{cc}\n", + "1 & 0 \\\\\n", + "0.5 & 1 \\end{array} \\right]\n", + "\\left(\\left[ \\begin{array}{cc}\n", + "2 & 1 \\\\\n", + "0 & 2.5 \\end{array} \\right]\n", + "\\left[\\begin{array}{c} \n", + "x_{1} \\\\ \n", + "x_{2} \\end{array}\\right]-\n", + "\\left[\\begin{array}{c} \n", + "1 \\\\\n", + "0.5\\end{array}\\right]\\right)=0$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "L =\n", + "\n", + " 1.00000 0.00000\n", + " 0.50000 1.00000\n", + "\n", + "U =\n", + "\n", + " 2.00000 1.00000\n", + " 0.00000 2.50000\n", + "\n", + "ans =\n", + "\n", + " 2 1\n", + " 1 3\n", + "\n", + "d =\n", + "\n", + " 1.00000\n", + " 0.50000\n", + "\n", + "y =\n", + "\n", + " 1\n", + " 1\n", + "\n" + ] + } + ], + "source": [ + "A=[2,1;1,3]\n", + "L=[1,0;0.5,1]\n", + "U=[2,1;0,2.5]\n", + "L*U\n", + "\n", + "d=[1;0.5]\n", + "y=L*d" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 5.0000\n", + "ans = 1\n", + "ans = 5\n", + "ans = 5\n", + "ans = 5.0000\n" + ] + } + ], + "source": [ + "% what is the determinant of L, U and A?\n", + "\n", + "det(A)\n", + "det(L)\n", + "det(U)\n", + "det(L)*det(U)\n", + "det(L*U)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pivoting for LU factorization\n", + "\n", + "LU factorization uses the same method as Gauss elimination so it is also necessary to perform partial pivoting when creating the lower and upper triangular matrices. \n", + "\n", + "Matlab and Octave use pivoting in the command \n", + "\n", + "`[L,U,P]=lu(A)`\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n", + "\n", + " -- Built-in Function: [L, U] = lu (A)\n", + " -- Built-in Function: [L, U, P] = lu (A)\n", + " -- Built-in Function: [L, U, P, Q] = lu (S)\n", + " -- Built-in Function: [L, U, P, Q, R] = lu (S)\n", + " -- Built-in Function: [...] = lu (S, THRES)\n", + " -- Built-in Function: Y = lu (...)\n", + " -- Built-in Function: [...] = lu (..., \"vector\")\n", + " Compute the LU decomposition of A.\n", + "\n", + " If A is full subroutines from LAPACK are used and if A is sparse\n", + " then UMFPACK is used.\n", + "\n", + " The result is returned in a permuted form, according to the\n", + " optional return value P. For example, given the matrix 'a = [1, 2;\n", + " 3, 4]',\n", + "\n", + " [l, u, p] = lu (A)\n", + "\n", + " returns\n", + "\n", + " l =\n", + "\n", + " 1.00000 0.00000\n", + " 0.33333 1.00000\n", + "\n", + " u =\n", + "\n", + " 3.00000 4.00000\n", + " 0.00000 0.66667\n", + "\n", + " p =\n", + "\n", + " 0 1\n", + " 1 0\n", + "\n", + " The matrix is not required to be square.\n", + "\n", + " When called with two or three output arguments and a spare input\n", + " matrix, 'lu' does not attempt to perform sparsity preserving column\n", + " permutations. Called with a fourth output argument, the sparsity\n", + " preserving column transformation Q is returned, such that 'P * A *\n", + " Q = L * U'.\n", + "\n", + " Called with a fifth output argument and a sparse input matrix, 'lu'\n", + " attempts to use a scaling factor R on the input matrix such that 'P\n", + " * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n", + " stable factorization.\n", + "\n", + " An additional input argument THRES, that defines the pivoting\n", + " threshold can be given. THRES can be a scalar, in which case it\n", + " defines the UMFPACK pivoting tolerance for both symmetric and\n", + " unsymmetric cases. If THRES is a 2-element vector, then the first\n", + " element defines the pivoting tolerance for the unsymmetric UMFPACK\n", + " pivoting strategy and the second for the symmetric strategy. By\n", + " default, the values defined by 'spparms' are used ([0.1, 0.001]).\n", + "\n", + " Given the string argument \"vector\", 'lu' returns the values of P\n", + " and Q as vector values, such that for full matrix, 'A (P,:) = L *\n", + " U', and 'R(P,:) * A (:, Q) = L * U'.\n", + "\n", + " With two output arguments, returns the permuted forms of the upper\n", + " and lower triangular matrices, such that 'A = L * U'. With one\n", + " output argument Y, then the matrix returned by the LAPACK routines\n", + " is returned. If the input matrix is sparse then the matrix L is\n", + " embedded into U to give a return value similar to the full case.\n", + " For both full and sparse matrices, 'lu' loses the permutation\n", + " information.\n", + "\n", + " See also: luupdate, ilu, chol, hess, qr, qz, schur, svd.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help lu" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tLU decomp\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tLU decomp\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tOctave \\\\\n", + "\n", + "\t\n", + "\t\tOctave \\\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% time LU solution vs backslash\n", + "t_lu=zeros(100,1);\n", + "t_bs=zeros(100,1);\n", + "for N=1:100\n", + " A=rand(N,N);\n", + " y=rand(N,1);\n", + " [L,U,P]=lu(A);\n", + "\n", + " tic; d=inv(L)*y; x=inv(U)*d; t_lu(N)=toc;\n", + "\n", + " tic; x=inv(A)*y; t_bs(N)=toc;\n", + "end\n", + "plot([1:100],t_lu,[1:100],t_bs) \n", + "legend('LU decomp','Octave \\\\')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n", + "\n", + "![Springs-masses](../lecture_09/mass_springs.svg)\n", + "\n", + "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n", + "\n", + "$m_{1}g+k(x_{2}-x_{1})-kx_{1}=0$\n", + "\n", + "$m_{2}g+k(x_{3}-x_{2})-k(x_{2}-x_{1})=0$\n", + "\n", + "$m_{3}g+k(x_{4}-x_{3})-k(x_{3}-x_{2})=0$\n", + "\n", + "$m_{4}g-k(x_{4}-x_{3})=0$\n", + "\n", + "in matrix form:\n", + "\n", + "$\\left[ \\begin{array}{cccc}\n", + "2k & -k & 0 & 0 \\\\\n", + "-k & 2k & -k & 0 \\\\\n", + "0 & -k & 2k & -k \\\\\n", + "0 & 0 & -k & k \\end{array} \\right]\n", + "\\left[ \\begin{array}{c}\n", + "x_{1} \\\\\n", + "x_{2} \\\\\n", + "x_{3} \\\\\n", + "x_{4} \\end{array} \\right]=\n", + "\\left[ \\begin{array}{c}\n", + "m_{1}g \\\\\n", + "m_{2}g \\\\\n", + "m_{3}g \\\\\n", + "m_{4}g \\end{array} \\right]$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 20 -10 0 0\n", + " -10 20 -10 0\n", + " 0 -10 20 -10\n", + " 0 0 -10 10\n", + "\n", + "y =\n", + "\n", + " 9.8100\n", + " 19.6200\n", + " 29.4300\n", + " 39.2400\n", + "\n" + ] + } + ], + "source": [ + "k=10; % N/m\n", + "m1=1; % kg\n", + "m2=2;\n", + "m3=3;\n", + "m4=4;\n", + "g=9.81; % m/s^2\n", + "K=[2*k -k 0 0; -k 2*k -k 0; 0 -k 2*k -k; 0 0 -k k]\n", + "y=[m1*g;m2*g;m3*g;m4*g]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This matrix, K, is symmetric. \n", + "\n", + "`K(i,j)==K(j,i)`\n", + "\n", + "Now we can use,\n", + "\n", + "## Cholesky Factorization\n", + "\n", + "We can decompose the matrix, K into two matrices, $U$ and $U^{T}$, where\n", + "\n", + "$K=U^{T}U$\n", + "\n", + "each of the components of U can be calculated with the following equations:\n", + "\n", + "$u_{ii}=\\sqrt{a_{ii}-\\sum_{k=1}^{i-1}u_{ki}^{2}}$\n", + "\n", + "$u_{ij}=\\frac{a_{ij}-\\sum_{k=1}^{i-1}u_{ki}u_{kj}}{u_{ii}}$\n", + "\n", + "so for K" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "K =\n", + "\n", + " 20 -10 0 0\n", + " -10 20 -10 0\n", + " 0 -10 20 -10\n", + " 0 0 -10 10\n", + "\n" + ] + } + ], + "source": [ + "K" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "u11 = 4.4721\n", + "u12 = -2.2361\n", + "u13 = 0\n", + "u14 = 0\n", + "u22 = 3.8730\n", + "u23 = -2.5820\n", + "u24 = 0\n", + "u33 = 3.6515\n", + "u34 = -2.7386\n", + "u44 = 1.5811\n", + "U =\n", + "\n", + " 4.47214 -2.23607 0.00000 0.00000\n", + " 0.00000 3.87298 -2.58199 0.00000\n", + " 0.00000 0.00000 3.65148 -2.73861\n", + " 0.00000 0.00000 0.00000 1.58114\n", + "\n" + ] + } + ], + "source": [ + "u11=sqrt(K(1,1))\n", + "u12=(K(1,2))/u11\n", + "u13=(K(1,3))/u11\n", + "u14=(K(1,4))/u11\n", + "u22=sqrt(K(2,2)-u12^2)\n", + "u23=(K(2,3)-u12*u13)/u22\n", + "u24=(K(2,4)-u12*u14)/u22\n", + "u33=sqrt(K(3,3)-u13^2-u23^2)\n", + "u34=(K(3,4)-u13*u14-u23*u24)/u33\n", + "u44=sqrt(K(4,4)-u14^2-u24^2-u34^2)\n", + "U=[u11,u12,u13,u14;0,u22,u23,u24;0,0,u33,u34;0,0,0,u44]" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans =\n", + "\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + " 1 1 1 1\n", + "\n" + ] + } + ], + "source": [ + "(U'*U)'==U'*U" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "average time spent for Cholesky factored solution = 1.465964e-05+/-9.806001e-06\n", + "average time spent for backslash solution = 1.555967e-05+/-1.048561e-05\n" + ] + } + ], + "source": [ + "% time solution for Cholesky vs backslash\n", + "t_chol=zeros(1000,1);\n", + "t_bs=zeros(1000,1);\n", + "for i=1:1000\n", + " tic; d=U'*y; x=U\\d; t_chol(i)=toc;\n", + " tic; x=K\\y; t_bs(i)=toc;\n", + "end\n", + "fprintf('average time spent for Cholesky factored solution = %e+/-%e',mean(t_chol),std(t_chol))\n", + "\n", + "fprintf('average time spent for backslash solution = %e+/-%e',mean(t_bs),std(t_bs))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}