Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
ME3255S2017/lecture_13/lecture_13.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
1401 lines (1401 sloc)
55.3 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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": [ | |
"<svg height=\"420px\" viewBox=\"0 0 560 420\" width=\"560px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n", | |
"\n", | |
"<title>Gnuplot</title>\n", | |
"<desc>Produced by GNUPLOT 5.0 patchlevel 3 </desc>\n", | |
"\n", | |
"<g id=\"gnuplot_canvas\">\n", | |
"\n", | |
"<rect fill=\"none\" height=\"420\" width=\"560\" x=\"0\" y=\"0\"/>\n", | |
"<defs>\n", | |
"\n", | |
"\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\"/>\n", | |
"\t<path d=\"M-1,0 h2 M0,-1 v2\" id=\"gpPt0\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<path d=\"M-1,-1 L1,1 M1,-1 L-1,1\" id=\"gpPt1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<path d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\" id=\"gpPt2\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<rect height=\"2\" id=\"gpPt3\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n", | |
"\t<rect fill=\"currentColor\" height=\"2\" id=\"gpPt4\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n", | |
"\t<circle cx=\"0\" cy=\"0\" id=\"gpPt5\" r=\"1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt6\" stroke=\"none\" xlink:href=\"#gpPt5\"/>\n", | |
"\t<path d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\" id=\"gpPt7\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt8\" stroke=\"none\" xlink:href=\"#gpPt7\"/>\n", | |
"\t<use id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\" xlink:href=\"#gpPt7\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt10\" stroke=\"none\" xlink:href=\"#gpPt9\"/>\n", | |
"\t<use id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\" xlink:href=\"#gpPt3\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt12\" stroke=\"none\" xlink:href=\"#gpPt11\"/>\n", | |
"\t<path d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\" id=\"gpPt13\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt14\" stroke=\"none\" xlink:href=\"#gpPt13\"/>\n", | |
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"textbox\" width=\"1\" x=\"0\" y=\"0\">\n", | |
"\t <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n", | |
"\t <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n", | |
"\t</filter>\n", | |
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"greybox\" width=\"1\" x=\"0\" y=\"0\">\n", | |
"\t <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n", | |
"\t <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n", | |
"\t</filter>\n", | |
"</defs>\n", | |
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"\t<g shape-rendering=\"crispEdges\" stroke=\"none\">\n", | |
"\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"70.5,384.0 534.9,384.0 534.9,16.8 70.5,16.8 \"/>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"rgb(255, 255, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,384.0 L83.0,384.0 M535.0,384.0 L522.5,384.0 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,390.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,292.2 L83.0,292.2 M535.0,292.2 L522.5,292.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,298.2)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.0005</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,200.3 L83.0,200.3 M535.0,200.3 L522.5,200.3 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,206.3)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.001</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,108.5 L83.0,108.5 M535.0,108.5 L522.5,108.5 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,114.5)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.0015</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,16.7 L83.0,16.7 M535.0,16.7 L522.5,16.7 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,22.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.002</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,384.0 L70.5,371.5 M70.5,16.7 L70.5,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(70.5,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M163.4,384.0 L163.4,371.5 M163.4,16.7 L163.4,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(163.4,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">20</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M256.3,384.0 L256.3,371.5 M256.3,16.7 L256.3,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(256.3,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">40</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M349.2,384.0 L349.2,371.5 M349.2,16.7 L349.2,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(349.2,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">60</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M442.1,384.0 L442.1,371.5 M442.1,16.7 L442.1,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(442.1,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">80</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M535.0,384.0 L535.0,371.5 M535.0,16.7 L535.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(535.0,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">100</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,16.7 L70.5,384.0 L535.0,384.0 L535.0,16.7 L70.5,16.7 Z \" stroke=\"black\"/></g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"\t<path d=\"M78.8,169.7 L78.8,25.7 L311.8,25.7 L311.8,169.7 L78.8,169.7 Z \" stroke=\"black\"/></g>\n", | |
"\t<g id=\"gnuplot_plot_1a\"><title>inversion</title>\n", | |
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(311.8,55.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">inversion</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<path d=\"M90.0,49.7 L143.8,49.7 M75.1,367.8 L79.8,374.3 L84.4,375.7 L89.1,375.5 L93.7,375.5 L98.4,375.2 L103.0,374.6 L107.7,374.5 L112.3,373.4 L117.0,372.8 L121.6,368.2 L126.2,372.4 L130.9,372.0 L135.5,371.5 L140.2,371.1 L144.8,369.7 L149.5,373.5 L154.1,250.8 L158.8,365.1 L163.4,364.9 L168.0,364.2 L172.7,362.5 L177.3,362.2 L182.0,361.1 L186.6,359.9 L191.3,360.1 L195.9,357.4 L200.6,358.3 L205.2,356.3 L209.9,362.5 L213.7,16.7 M215.3,16.7 L219.1,360.7 L223.8,360.5 L228.4,356.3 L233.1,354.6 L237.7,356.8 L242.4,354.4 L247.0,360.1 L251.7,359.6 L256.3,359.0 L260.9,358.5 L265.6,356.1 L270.2,356.5 L274.9,354.8 L279.5,354.1 L284.2,351.9 L288.8,351.3 L293.5,352.0 L298.1,350.0 L302.8,348.4 L307.4,348.0 L312.0,354.4 L316.7,351.5 L321.3,352.4 L326.0,351.3 L330.6,351.7 L335.3,346.3 L339.9,345.4 L344.6,347.3 L349.2,347.4 L353.8,345.2 L358.5,342.5 L363.1,344.3 L367.8,347.1 L372.4,338.9 L377.1,340.1 L381.7,340.1 L386.4,329.3 L391.0,337.0 L395.7,319.7 L400.3,338.6 L404.9,340.0 L409.6,336.3 L414.2,335.3 L418.9,333.7 L423.5,333.3 L428.2,330.2 L432.8,328.7 L437.5,332.8 L442.1,327.4 L446.7,325.6 L451.4,326.5 L456.0,325.2 L460.7,324.5 L465.3,320.3 L470.0,325.4 L474.6,319.4 L479.3,322.1 L483.9,319.7 L488.6,321.2 L493.2,315.7 L497.8,316.6 L502.5,299.5 L507.1,317.5 L511.8,310.4 L516.4,318.2 L521.1,314.0 L525.7,316.6 L530.4,311.8 L535.0,314.0 \" stroke=\"rgb( 0, 0, 255)\"/></g>\n", | |
"\t</g>\n", | |
"\t<g id=\"gnuplot_plot_2a\"><title>backslash</title>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(311.8,103.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">backslash</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<path d=\"M90.0,97.7 L143.8,97.7 M75.1,381.3 L79.8,377.7 L84.4,378.8 L89.1,378.7 L93.7,378.7 L98.4,377.2 L103.0,378.0 L107.7,378.0 L112.3,377.2 L117.0,376.6 L121.6,376.3 L126.2,376.1 L130.9,375.8 L135.5,375.9 L140.2,375.5 L144.8,377.6 L149.5,377.0 L154.1,370.4 L158.8,371.0 L163.4,371.0 L168.0,371.0 L172.7,369.3 L177.3,370.0 L182.0,369.3 L186.6,368.6 L191.3,368.0 L195.9,367.3 L200.6,366.6 L205.2,366.7 L209.9,370.3 L214.5,367.7 L219.1,369.5 L223.8,340.8 L228.4,367.3 L233.1,367.6 L237.7,367.1 L242.4,368.9 L247.0,369.3 L251.7,368.6 L256.3,368.8 L260.9,368.4 L265.6,367.7 L270.2,367.1 L274.9,365.3 L279.5,366.6 L284.2,364.2 L288.8,365.1 L293.5,365.6 L298.1,364.9 L302.8,363.1 L307.4,366.0 L312.0,366.2 L316.7,363.8 L321.3,365.6 L326.0,365.1 L330.6,365.3 L335.3,362.5 L339.9,362.3 L344.6,363.4 L349.2,363.1 L353.8,362.9 L358.5,364.0 L363.1,361.1 L367.8,363.8 L372.4,361.4 L377.1,352.7 L381.7,357.9 L386.4,359.7 L391.0,355.2 L395.7,358.3 L400.3,357.0 L404.9,358.5 L409.6,357.6 L414.2,356.8 L418.9,354.6 L423.5,355.1 L428.2,355.9 L432.8,353.1 L437.5,351.5 L442.1,351.9 L446.7,351.7 L451.4,351.1 L456.0,350.9 L460.7,350.5 L465.3,347.7 L470.0,350.4 L474.6,349.5 L479.3,349.3 L483.9,348.9 L488.6,348.4 L493.2,342.5 L497.8,346.5 L502.5,347.8 L507.1,347.8 L511.8,347.1 L516.4,347.6 L521.1,347.3 L525.7,346.2 L530.4,345.6 L535.0,346.0 \" stroke=\"rgb( 0, 128, 0)\"/></g>\n", | |
"\t</g>\n", | |
"\t<g id=\"gnuplot_plot_3a\"><title>multiplication</title>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(311.8,151.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">multiplication</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<path d=\"M90.0,145.7 L143.8,145.7 M75.1,382.0 L79.8,381.4 L84.4,381.8 L89.1,381.8 L93.7,381.8 L98.4,381.8 L103.0,381.8 L107.7,381.8 L112.3,381.8 L117.0,381.9 L121.6,381.8 L126.2,381.8 L130.9,382.0 L135.5,381.8 L140.2,381.8 L144.8,382.6 L149.5,382.3 L154.1,381.8 L158.8,381.8 L163.4,382.0 L168.0,381.8 L172.7,381.8 L177.3,381.8 L182.0,381.8 L186.6,381.8 L191.3,381.8 L195.9,381.6 L200.6,381.6 L205.2,381.8 L209.9,382.2 L214.5,381.9 L219.1,382.2 L223.8,382.0 L228.4,382.2 L233.1,382.2 L237.7,382.0 L242.4,382.6 L247.0,382.3 L251.7,382.6 L256.3,382.6 L260.9,382.4 L265.6,382.5 L270.2,382.5 L274.9,382.3 L279.5,382.3 L284.2,382.3 L288.8,382.5 L293.5,382.4 L298.1,382.3 L302.8,379.9 L307.4,382.7 L312.0,382.5 L316.7,382.5 L321.3,382.7 L326.0,382.7 L330.6,382.6 L335.3,382.6 L339.9,382.5 L344.6,382.6 L349.2,382.5 L353.8,382.6 L358.5,382.7 L363.1,382.6 L367.8,382.7 L372.4,382.5 L377.1,382.6 L381.7,382.5 L386.4,382.5 L391.0,382.6 L395.7,382.5 L400.3,382.5 L404.9,382.5 L409.6,382.6 L414.2,382.5 L418.9,382.5 L423.5,382.6 L428.2,382.6 L432.8,382.6 L437.5,382.3 L442.1,382.5 L446.7,382.6 L451.4,382.5 L456.0,382.5 L460.7,382.5 L465.3,382.4 L470.0,382.3 L474.6,382.3 L479.3,382.6 L483.9,382.4 L488.6,382.3 L493.2,382.3 L497.8,382.6 L502.5,382.6 L507.1,382.3 L511.8,382.3 L516.4,382.3 L521.1,382.6 L525.7,382.3 L530.4,382.3 L535.0,382.4 \" stroke=\"rgb(255, 0, 0)\"/></g>\n", | |
"\t</g>\n", | |
"<g color=\"white\" fill=\"none\" stroke=\"rgb(255, 0, 0)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"</g>\n", | |
"</svg>" | |
], | |
"text/plain": [ | |
"<IPython.core.display.SVG object>" | |
] | |
}, | |
"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 <topic>' 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": [ | |
"<svg height=\"420px\" viewBox=\"0 0 560 420\" width=\"560px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n", | |
"\n", | |
"<title>Gnuplot</title>\n", | |
"<desc>Produced by GNUPLOT 5.0 patchlevel 3 </desc>\n", | |
"\n", | |
"<g id=\"gnuplot_canvas\">\n", | |
"\n", | |
"<rect fill=\"none\" height=\"420\" width=\"560\" x=\"0\" y=\"0\"/>\n", | |
"<defs>\n", | |
"\n", | |
"\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\"/>\n", | |
"\t<path d=\"M-1,0 h2 M0,-1 v2\" id=\"gpPt0\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<path d=\"M-1,-1 L1,1 M1,-1 L-1,1\" id=\"gpPt1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<path d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\" id=\"gpPt2\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<rect height=\"2\" id=\"gpPt3\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n", | |
"\t<rect fill=\"currentColor\" height=\"2\" id=\"gpPt4\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n", | |
"\t<circle cx=\"0\" cy=\"0\" id=\"gpPt5\" r=\"1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt6\" stroke=\"none\" xlink:href=\"#gpPt5\"/>\n", | |
"\t<path d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\" id=\"gpPt7\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt8\" stroke=\"none\" xlink:href=\"#gpPt7\"/>\n", | |
"\t<use id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\" xlink:href=\"#gpPt7\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt10\" stroke=\"none\" xlink:href=\"#gpPt9\"/>\n", | |
"\t<use id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\" xlink:href=\"#gpPt3\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt12\" stroke=\"none\" xlink:href=\"#gpPt11\"/>\n", | |
"\t<path d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\" id=\"gpPt13\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n", | |
"\t<use fill=\"currentColor\" id=\"gpPt14\" stroke=\"none\" xlink:href=\"#gpPt13\"/>\n", | |
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"textbox\" width=\"1\" x=\"0\" y=\"0\">\n", | |
"\t <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n", | |
"\t <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n", | |
"\t</filter>\n", | |
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"greybox\" width=\"1\" x=\"0\" y=\"0\">\n", | |
"\t <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n", | |
"\t <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n", | |
"\t</filter>\n", | |
"</defs>\n", | |
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"\t<g shape-rendering=\"crispEdges\" stroke=\"none\">\n", | |
"\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"70.5,384.0 534.9,384.0 534.9,16.8 70.5,16.8 \"/>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"rgb(255, 255, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,384.0 L83.0,384.0 M535.0,384.0 L522.5,384.0 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,390.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,292.2 L83.0,292.2 M535.0,292.2 L522.5,292.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,298.2)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.0005</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,200.3 L83.0,200.3 M535.0,200.3 L522.5,200.3 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,206.3)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.001</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,108.5 L83.0,108.5 M535.0,108.5 L522.5,108.5 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,114.5)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.0015</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,16.7 L83.0,16.7 M535.0,16.7 L522.5,16.7 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(62.2,22.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0.002</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,384.0 L70.5,371.5 M70.5,16.7 L70.5,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(70.5,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M163.4,384.0 L163.4,371.5 M163.4,16.7 L163.4,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(163.4,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">20</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M256.3,384.0 L256.3,371.5 M256.3,16.7 L256.3,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(256.3,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">40</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M349.2,384.0 L349.2,371.5 M349.2,16.7 L349.2,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(349.2,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">60</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M442.1,384.0 L442.1,371.5 M442.1,16.7 L442.1,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(442.1,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">80</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M535.0,384.0 L535.0,371.5 M535.0,16.7 L535.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(535.0,408.0)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">100</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"\t<path d=\"M70.5,16.7 L70.5,384.0 L535.0,384.0 L535.0,16.7 L70.5,16.7 Z \" stroke=\"black\"/></g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n", | |
"\t<path d=\"M349.7,121.7 L349.7,25.7 L526.7,25.7 L526.7,121.7 L349.7,121.7 Z \" stroke=\"black\"/></g>\n", | |
"\t<g id=\"gnuplot_plot_1a\"><title>LU decomp</title>\n", | |
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(526.7,55.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">LU decomp</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<path d=\"M360.9,49.7 L414.7,49.7 M75.1,357.4 L79.8,366.4 L84.4,370.3 L89.1,369.9 L93.7,363.8 L98.4,369.5 L103.0,369.5 L107.7,362.2 L112.3,368.4 L117.0,368.0 L121.6,360.8 L126.2,366.5 L130.9,366.6 L135.5,355.4 L140.2,355.5 L144.8,353.0 L149.5,346.7 L154.1,360.0 L158.8,362.3 L163.4,361.2 L168.0,362.3 L172.7,355.7 L177.3,360.3 L182.0,356.9 L186.6,355.1 L191.3,357.9 L195.9,351.3 L200.6,356.6 L205.2,355.0 L209.9,355.4 L214.5,362.2 L219.1,360.9 L223.8,360.8 L228.4,359.7 L233.1,359.0 L237.7,356.5 L242.4,358.3 L247.0,356.6 L251.7,360.7 L256.3,360.3 L260.9,360.0 L265.6,359.0 L270.2,359.0 L274.9,357.6 L279.5,357.4 L284.2,356.2 L288.8,354.2 L293.5,359.6 L298.1,358.5 L302.8,357.9 L307.4,356.5 L312.0,356.5 L316.7,355.9 L321.3,355.1 L326.0,353.9 L330.6,352.4 L335.3,357.0 L339.9,356.6 L344.6,355.9 L349.2,355.2 L353.8,354.1 L358.5,353.1 L363.1,352.6 L367.8,351.5 L372.4,335.7 L377.1,356.3 L381.7,356.6 L386.4,352.6 L391.0,355.5 L395.7,355.5 L400.3,354.6 L404.9,355.1 L409.6,224.4 L414.2,352.4 L418.9,338.8 L423.5,351.5 L428.2,350.8 L432.8,350.4 L437.5,349.7 L442.1,350.4 L446.7,348.4 L451.4,347.8 L456.0,346.2 L460.7,347.0 L465.3,340.7 L470.0,338.6 L474.6,346.2 L479.3,347.3 L483.9,344.2 L488.6,343.8 L493.2,342.5 L497.8,343.2 L502.5,341.2 L507.1,342.5 L511.8,338.8 L516.4,342.3 L521.1,337.9 L525.7,338.3 L530.4,337.8 L535.0,339.2 \" stroke=\"rgb( 0, 0, 255)\"/></g>\n", | |
"\t</g>\n", | |
"\t<g id=\"gnuplot_plot_2a\"><title>Octave \\\\</title>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(526.7,103.7)\">\n", | |
"\t\t<text><tspan font-family=\"{}\">Octave \\</tspan></text>\n", | |
"\t</g>\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n", | |
"\t<path d=\"M360.9,97.7 L414.7,97.7 M75.1,371.7 L79.8,368.6 L84.4,372.8 L89.1,372.4 L93.7,371.3 L98.4,371.7 L103.0,371.3 L107.7,369.9 L112.3,370.3 L117.0,370.2 L121.6,368.6 L126.2,368.9 L130.9,364.9 L135.5,357.0 L140.2,358.5 L144.8,354.4 L149.5,353.0 L154.1,358.7 L158.8,359.8 L163.4,360.1 L168.0,360.7 L172.7,355.2 L177.3,353.7 L182.0,355.4 L186.6,337.5 L191.3,353.7 L195.9,349.1 L200.6,354.4 L205.2,350.4 L209.9,329.8 L214.5,357.6 L219.1,353.7 L223.8,357.6 L228.4,350.4 L233.1,354.4 L237.7,353.9 L242.4,352.4 L247.0,140.3 L251.7,355.1 L256.3,356.8 L260.9,355.5 L265.6,352.6 L270.2,353.3 L274.9,353.9 L279.5,351.7 L284.2,350.8 L288.8,341.2 L293.5,355.2 L298.1,353.0 L302.8,352.6 L307.4,350.6 L312.0,351.3 L316.7,343.8 L321.3,348.0 L326.0,347.3 L330.6,345.4 L335.3,351.5 L339.9,338.1 L344.6,348.7 L349.2,349.1 L353.8,347.3 L358.5,346.9 L363.1,343.2 L367.8,343.6 L372.4,343.2 L377.1,335.7 L381.7,342.7 L386.4,344.5 L391.0,345.8 L395.7,345.6 L400.3,332.4 L404.9,340.7 L409.6,341.9 L414.2,262.6 L418.9,340.1 L423.5,331.8 L428.2,339.4 L432.8,338.3 L437.5,334.6 L442.1,335.5 L446.7,328.7 L451.4,333.3 L456.0,328.7 L460.7,331.3 L465.3,329.7 L470.0,327.8 L474.6,323.0 L479.3,324.8 L483.9,114.6 L488.6,315.3 L493.2,314.0 L497.8,301.7 L502.5,315.3 L507.1,319.2 L511.8,314.6 L516.4,315.1 L521.1,311.6 L525.7,313.1 L530.4,302.3 L535.0,305.2 \" stroke=\"rgb( 0, 128, 0)\"/></g>\n", | |
"\t</g>\n", | |
"<g color=\"white\" fill=\"none\" stroke=\"rgb( 0, 128, 0)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n", | |
"</g>\n", | |
"</g>\n", | |
"</svg>" | |
], | |
"text/plain": [ | |
"<IPython.core.display.SVG object>" | |
] | |
}, | |
"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 | |
} |