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.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/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..e69de29 100644 Binary files a/lecture_11/octave-workspace and b/lecture_11/octave-workspace differ diff --git a/lecture_13/.ipynb_checkpoints/lecture_13-checkpoint.ipynb b/lecture_13/.ipynb_checkpoints/lecture_13-checkpoint.ipynb new file mode 100644 index 0000000..6a5dc68 --- /dev/null +++ b/lecture_13/.ipynb_checkpoints/lecture_13-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 +} diff --git a/lecture_13/lecture_13.ipynb b/lecture_13/lecture_13.ipynb new file mode 100644 index 0000000..c687623 --- /dev/null +++ b/lecture_13/lecture_13.ipynb @@ -0,0 +1,1401 @@ +{ + "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 t" + ] + }, + { + "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_13/octave-workspace b/lecture_13/octave-workspace new file mode 100644 index 0000000..e69de29