Skip to content
Permalink
master
Switch branches/tags

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_12/lecture_12.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
1248 lines (1248 sloc) 39.8 KB
Raw Blame
Open in GitHub Desktop
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%plot --format svg"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"setdefaults"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A =\n",
"\n",
" 0.447394 0.357071 0.720915 0.499926\n",
" 0.648313 0.323276 0.521677 0.288345\n",
" 0.084982 0.581513 0.466420 0.142342\n",
" 0.576580 0.658089 0.916987 0.923165\n",
"\n",
"L =\n",
"\n",
" 1.00000 0.00000 0.00000 0.00000\n",
" 0.13108 1.00000 0.00000 0.00000\n",
" 0.69009 0.24851 1.00000 0.00000\n",
" 0.88935 0.68736 0.68488 1.00000\n",
"\n",
"U =\n",
"\n",
" 0.64831 0.32328 0.52168 0.28834\n",
" 0.00000 0.53914 0.39804 0.10455\n",
" 0.00000 0.00000 0.26199 0.27496\n",
" 0.00000 0.00000 0.00000 0.40655\n",
"\n",
"P =\n",
"\n",
"Permutation Matrix\n",
"\n",
" 0 1 0 0\n",
" 0 0 1 0\n",
" 1 0 0 0\n",
" 0 0 0 1\n",
"\n",
"ans = 1\n"
]
}
],
"source": [
"A=rand(4,4)\n",
"\n",
"[L,U,P]=lu(A)\n",
"\n",
"det(L)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 4 4\n",
"\n",
"ans = 23.586\n",
"ans = 35.826\n",
"ans = 14.869\n",
"C =\n",
"\n",
" 5.98549 4.28555 4.35707 4.31359\n",
" 0.00000 3.63950 1.35005 1.45342\n",
" 0.00000 0.00000 3.62851 1.50580\n",
" 0.00000 0.00000 0.00000 3.21911\n",
"\n"
]
}
],
"source": [
"A=rand(4,100)';\n",
"A=A'*A;\n",
"size(A)\n",
"min(min(A))\n",
"max(max(A))\n",
"cond(A)\n",
"C=chol(A)"
]
},
{
"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": 45,
"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_{n} \\\\\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": "markdown",
"metadata": {},
"source": [
"#### Note on solving for $A^{-1}$ column 1\n",
"\n",
"$Aa_1=I(:,1)$\n",
"\n",
"$LUa_1=I(:,1)$\n",
"\n",
"$(LUa_1-I(:,1))=0$\n",
"\n",
"$L(Ua_1-d_1)=0$\n",
"\n",
"$I(:,1)=Ld_1$"
]
},
{
"cell_type": "code",
"execution_count": 56,
"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": 57,
"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",
"\\left[\\begin{array}{ccc} \n",
"1 & 0 & 0 \\\\ \n",
"-1/2 & 1 & 0 \\\\\n",
"0 & -2/3 & 1 \\end{array} \\right]\\left[\\begin{array}{c} \n",
"d1(1) \\\\ \n",
"d1(2) \\\\ \n",
"d1(3)\\end{array} \\right]=\\left[\\begin{array}{c} \n",
"1 \\\\ \n",
"0 \\\\ \n",
"0 \\end{array} \\right]\n",
";~Ua_{1}=d_{1}$"
]
},
{
"cell_type": "code",
"execution_count": 58,
"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": 59,
"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": 60,
"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": 61,
"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": 62,
"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": 63,
"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": 69,
"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",
"I_app =\n",
"\n",
" 1.00000 0.00000 0.00000\n",
" 0.00000 1.00000 -0.00000\n",
" -0.00000 -0.00000 1.00000\n",
"\n",
"ans = -4.4409e-16\n",
"ans = 2.2204e-16\n",
"ans = 0.0039062\n"
]
}
],
"source": [
"invA=[a1,a2,a3]\n",
"I_app=A*invA\n",
"I_app(2,3)\n",
"eps\n",
"\n",
"2^-8"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the solution of $x$ to $Ax=y$ is $x=A^{-1}y$"
]
},
{
"cell_type": "code",
"execution_count": 70,
"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": 4,
"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,365.8 L79.8,371.5 L84.4,373.0 L89.1,373.5 L93.7,373.1 L98.4,372.6 L103.0,372.2 L107.7,371.9 L112.3,371.1 L117.0,371.3 L121.6,370.2 L126.2,370.6 L130.9,367.8 L135.5,368.4 L140.2,367.8 L144.8,366.9 L149.5,366.4 L154.1,332.1 L158.8,359.6 L163.4,361.2 L168.0,360.5 L172.7,360.3 L177.3,357.2 L182.0,357.6 L186.6,351.1 L191.3,355.9 L195.9,351.5 L200.6,353.7 L205.2,351.7 L209.9,351.3 L212.7,16.7 M216.4,16.7 L219.1,330.9 L223.8,327.4 L228.4,324.3 L233.1,326.2 L237.7,323.8 L242.4,322.5 L247.0,319.7 L251.7,316.8 L256.3,317.0 L260.9,312.9 L265.6,311.3 L270.2,304.0 L274.9,308.9 L279.5,306.3 L284.2,329.1 L288.8,327.3 L293.5,327.4 L298.1,323.6 L302.8,336.6 L307.4,336.1 L312.0,336.6 L316.7,333.9 L321.3,331.9 L326.0,330.4 L330.6,331.5 L335.3,326.5 L339.9,337.7 L344.6,335.9 L349.2,336.6 L353.8,333.3 L358.5,327.6 L363.1,331.2 L367.8,328.9 L372.4,327.4 L377.1,335.1 L381.7,333.5 L386.4,336.3 L391.0,332.8 L395.7,332.2 L400.3,331.1 L404.9,331.3 L409.6,328.0 L414.2,336.4 L418.9,331.6 L423.5,334.8 L428.2,333.9 L432.8,329.7 L437.5,329.6 L442.1,331.3 L446.7,325.6 L451.4,332.9 L456.0,332.4 L460.7,332.4 L465.3,328.0 L470.0,329.8 L474.6,322.5 L479.3,326.9 L483.9,325.8 L488.6,322.3 L493.2,322.5 L497.8,323.4 L502.5,322.1 L507.1,321.9 L511.8,319.4 L516.4,322.7 L521.1,318.5 L525.7,315.9 L530.4,315.0 L535.0,312.2 \" 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,380.8 L79.8,375.5 L84.4,377.0 L89.1,377.2 L93.7,377.2 L98.4,376.9 L103.0,376.5 L107.7,376.5 L112.3,375.9 L117.0,376.1 L121.6,375.2 L126.2,374.8 L130.9,373.0 L135.5,373.7 L140.2,373.2 L144.8,372.4 L149.5,372.2 L154.1,368.0 L158.8,368.6 L163.4,368.4 L168.0,367.7 L172.7,367.6 L177.3,365.4 L182.0,366.4 L186.6,363.1 L191.3,364.9 L195.9,362.9 L200.6,363.2 L205.2,362.5 L209.9,362.2 L214.5,353.0 L219.1,351.7 L223.8,347.4 L228.4,347.3 L233.1,347.3 L237.7,345.8 L242.4,344.7 L247.0,344.1 L251.7,342.3 L256.3,342.1 L260.9,340.7 L265.6,339.2 L270.2,333.7 L274.9,336.8 L279.5,351.7 L284.2,350.9 L288.8,350.2 L292.9,16.7 M294.0,16.7 L298.1,357.4 L302.8,355.5 L307.4,355.7 L312.0,355.5 L316.7,353.7 L321.3,354.2 L326.0,352.4 L330.6,353.4 L335.3,357.9 L339.9,357.6 L344.6,356.6 L349.2,356.9 L353.8,355.9 L358.5,353.2 L363.1,355.0 L367.8,354.8 L372.4,358.5 L377.1,355.7 L381.7,353.0 L386.4,355.5 L391.0,354.6 L395.7,353.7 L400.3,353.4 L404.9,353.7 L409.6,356.8 L414.2,356.5 L418.9,355.7 L423.5,355.7 L428.2,355.0 L432.8,352.3 L437.5,353.7 L442.1,352.4 L446.7,355.3 L451.4,355.2 L456.0,354.1 L460.7,353.7 L465.3,350.0 L470.0,353.2 L474.6,352.0 L479.3,350.0 L483.9,351.7 L488.6,348.9 L493.2,345.8 L497.8,350.9 L502.5,349.7 L507.1,349.5 L511.8,348.9 L516.4,349.5 L521.1,348.8 L525.7,346.3 L530.4,346.7 L535.0,343.8 \" 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,381.6 L79.8,381.1 L84.4,381.2 L89.1,381.2 L93.7,381.4 L98.4,381.2 L103.0,381.2 L107.7,381.4 L112.3,381.4 L117.0,381.2 L121.6,381.2 L126.2,381.4 L130.9,381.4 L135.5,381.2 L140.2,381.2 L144.8,381.2 L149.5,381.2 L154.1,380.9 L158.8,381.2 L163.4,381.2 L168.0,381.1 L172.7,381.2 L177.3,381.2 L182.0,381.2 L186.6,381.1 L191.3,381.1 L195.9,381.2 L200.6,378.1 L205.2,381.2 L209.9,381.1 L214.5,379.4 L219.1,380.0 L223.8,379.9 L228.4,380.1 L233.1,380.0 L237.7,380.1 L242.4,380.0 L247.0,379.9 L251.7,380.0 L256.3,379.8 L260.9,379.8 L265.6,380.0 L270.2,379.8 L274.9,379.8 L279.5,381.2 L284.2,381.2 L288.8,381.1 L293.5,381.1 L298.1,382.0 L302.8,381.8 L307.4,381.8 L312.0,381.8 L316.7,381.8 L321.3,381.8 L326.0,381.8 L330.6,381.8 L335.3,382.2 L339.9,382.2 L344.6,382.2 L349.2,382.2 L353.8,382.2 L358.5,381.9 L363.1,382.2 L367.8,382.2 L372.4,382.3 L377.1,382.5 L381.7,382.3 L386.4,382.4 L391.0,382.4 L395.7,382.3 L400.3,382.3 L404.9,381.6 L409.6,382.3 L414.2,382.5 L418.9,382.6 L423.5,382.6 L428.2,382.6 L432.8,382.5 L437.5,382.5 L442.1,382.5 L446.7,382.7 L451.4,382.7 L456.0,382.7 L460.7,382.5 L465.3,382.5 L470.0,382.6 L474.6,382.5 L479.3,382.7 L483.9,382.7 L488.6,382.6 L493.2,382.7 L497.8,382.6 L502.5,382.5 L507.1,382.6 L511.8,382.5 L516.4,382.4 L521.1,382.3 L525.7,382.5 L530.4,382.3 L535.0,382.3 \" 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 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": 72,
"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": 75,
"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": 74,
"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
}