diff --git a/lecture_13/lecture_13.ipynb b/lecture_12/.ipynb_checkpoints/lecture_12-checkpoint.ipynb similarity index 96% rename from lecture_13/lecture_13.ipynb rename to lecture_12/.ipynb_checkpoints/lecture_12-checkpoint.ipynb index c687623..1332b29 100644 --- a/lecture_13/lecture_13.ipynb +++ b/lecture_12/.ipynb_checkpoints/lecture_12-checkpoint.ipynb @@ -768,7 +768,89 @@ "### *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" + "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=" ] }, { 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_13-checkpoint.ipynb b/lecture_13/.ipynb_checkpoints/lecture_12-checkpoint.ipynb similarity index 100% rename from lecture_13/.ipynb_checkpoints/lecture_13-checkpoint.ipynb rename to lecture_13/.ipynb_checkpoints/lecture_12-checkpoint.ipynb diff --git a/lecture_13/octave-workspace b/lecture_13/octave-workspace deleted file mode 100644 index e69de29..0000000