diff --git a/09_Linear-Algebra/.ipynb_checkpoints/09_Linear-Algebra-checkpoint.ipynb b/09_Linear-Algebra/.ipynb_checkpoints/09_Linear-Algebra-checkpoint.ipynb
index 6bc2ad3..cba99b5 100644
--- a/09_Linear-Algebra/.ipynb_checkpoints/09_Linear-Algebra-checkpoint.ipynb
+++ b/09_Linear-Algebra/.ipynb_checkpoints/09_Linear-Algebra-checkpoint.ipynb
@@ -205,7 +205,7 @@
},
{
"cell_type": "code",
- "execution_count": 4,
+ "execution_count": 7,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -227,18 +227,25 @@
" 1 2 3 4\n",
" 5 6 7 8\n",
" 9 10 11 12\n",
+ "\n",
+ "C =\n",
+ "\n",
+ " 20 28 36 44\n",
+ " 38 44 50 56\n",
"\n"
]
}
],
"source": [
"A=[5,3,0;1,2,3] \n",
- "B=[1,2,3,4;5,6,7,8;9,10,11,12]"
+ "B=[1,2,3,4;5,6,7,8;9,10,11,12]\n",
+ "C=A*B;\n",
+ "C"
]
},
{
"cell_type": "code",
- "execution_count": 7,
+ "execution_count": 6,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -372,10 +379,19 @@
"\n",
"Consider 4 masses connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
"\n",
- "\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",
"\n",
+ "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"$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",
@@ -414,7 +430,7 @@
},
{
"cell_type": "code",
- "execution_count": 11,
+ "execution_count": 10,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -426,20 +442,6 @@
"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",
"x =\n",
"\n",
" 9.8100\n",
@@ -457,17 +459,21 @@
"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]\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];\n",
"\n",
"x=K\\y"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
- "## Matrix Operations\n",
+ "## Special Matrices\n",
"\n",
"Identity matrix `eye(M,N)` **Assume M=N unless specfied**\n",
"\n",
@@ -487,8 +493,13 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
"source": [
+ "## Matrix operations\n",
"### Transpose\n",
"\n",
"The transpose of a matrix changes the rows -> columns and columns-> rows\n",
@@ -514,9 +525,12 @@
},
{
"cell_type": "code",
- "execution_count": 10,
+ "execution_count": 13,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
},
"outputs": [
{
@@ -602,7 +616,7 @@
},
{
"cell_type": "code",
- "execution_count": 11,
+ "execution_count": 15,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -619,15 +633,16 @@
"Diagonal Matrix\n",
"\n",
" 1 0 0\n",
- " 0 1 0\n",
+ " 0 4 0\n",
" 0 0 1\n",
"\n",
- "ans = 3\n"
+ "ans = 6\n"
]
}
],
"source": [
- "id_m=eye(3)\n",
+ "id_m=eye(3);\n",
+ "id_m(2,2)=4\n",
"trace(id_m)"
]
},
@@ -658,7 +673,7 @@
},
{
"cell_type": "code",
- "execution_count": 20,
+ "execution_count": 16,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -672,9 +687,9 @@
"text": [
"A =\n",
"\n",
- " 0.5762106 0.3533174 0.7172134\n",
- " 0.7061664 0.4863733 0.9423436\n",
- " 0.4255961 0.0016613 0.3561407\n",
+ " 0.71624 0.28713 0.91248\n",
+ " 0.43423 0.50924 0.89731\n",
+ " 0.81776 0.15682 0.48825\n",
"\n"
]
}
@@ -685,7 +700,7 @@
},
{
"cell_type": "code",
- "execution_count": 22,
+ "execution_count": 17,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -699,9 +714,9 @@
"text": [
"Ainv =\n",
"\n",
- " 41.5613 -30.1783 -3.8467\n",
- " 36.2130 -24.2201 -8.8415\n",
- " -49.8356 36.1767 7.4460\n",
+ " -1.189272 -0.032028 2.281478\n",
+ " -5.750167 4.369462 2.716129\n",
+ " 3.838847 -1.349811 -2.645516\n",
"\n"
]
}
@@ -712,7 +727,7 @@
},
{
"cell_type": "code",
- "execution_count": 26,
+ "execution_count": 18,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -726,9 +741,9 @@
"text": [
"B =\n",
"\n",
- " 0.524529 0.470856 0.708116\n",
- " 0.084491 0.730986 0.528292\n",
- " 0.303545 0.782522 0.389282\n",
+ " 0.88372 0.43914 0.13476\n",
+ " 0.43439 0.34519 0.30337\n",
+ " 0.84604 0.86479 0.21558\n",
"\n"
]
}
@@ -739,7 +754,7 @@
},
{
"cell_type": "code",
- "execution_count": 28,
+ "execution_count": 20,
"metadata": {
"collapsed": false,
"slideshow": {
@@ -753,34 +768,22 @@
"text": [
"ans =\n",
"\n",
- " -182.185 125.738 40.598\n",
- " -133.512 97.116 17.079\n",
- " 282.422 -200.333 -46.861\n",
- "\n",
- "ans =\n",
- "\n",
- " -182.185 125.738 40.598\n",
- " -133.512 97.116 17.079\n",
- " 282.422 -200.333 -46.861\n",
- "\n",
- "ans =\n",
- "\n",
- " 41.5613 36.2130 -49.8356\n",
- " -30.1783 -24.2201 36.1767\n",
- " -3.8467 -8.8415 7.4460\n",
+ " -1.189272 -5.750167 3.838847\n",
+ " -0.032028 4.369462 -1.349811\n",
+ " 2.281478 2.716129 -2.645516\n",
"\n",
"ans =\n",
"\n",
- " 41.5613 36.2130 -49.8356\n",
- " -30.1783 -24.2201 36.1767\n",
- " -3.8467 -8.8415 7.4460\n",
+ " -1.189272 -5.750167 3.838847\n",
+ " -0.032028 4.369462 -1.349811\n",
+ " 2.281478 2.716129 -2.645516\n",
"\n"
]
}
],
"source": [
- "inv(A*B)\n",
- "inv(B)*inv(A)\n",
+ "%inv(A*B)\n",
+ "%inv(B)*inv(A)\n",
"\n",
"inv(A')\n",
"\n",
diff --git a/09_Linear-Algebra/09_Linear-Algebra.ipynb b/09_Linear-Algebra/09_Linear-Algebra.ipynb
index 71dc415..cba99b5 100644
--- a/09_Linear-Algebra/09_Linear-Algebra.ipynb
+++ b/09_Linear-Algebra/09_Linear-Algebra.ipynb
@@ -379,7 +379,7 @@
"\n",
"Consider 4 masses connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
"\n",
- "\n",
+ "\n",
"\n",
"The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:"
]
diff --git a/09_Linear-Algebra/octave-workspace b/09_Linear-Algebra/octave-workspace
index e69de29..8c437bb 100644
Binary files a/09_Linear-Algebra/octave-workspace and b/09_Linear-Algebra/octave-workspace differ
diff --git a/10_Gauss_elimination/.ipynb_checkpoints/10_Gauss_elimination-checkpoint.ipynb b/10_Gauss_elimination/.ipynb_checkpoints/10_Gauss_elimination-checkpoint.ipynb
new file mode 100644
index 0000000..1f6b0bb
--- /dev/null
+++ b/10_Gauss_elimination/.ipynb_checkpoints/10_Gauss_elimination-checkpoint.ipynb
@@ -0,0 +1,1802 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Gauss Elimination\n",
+ "### Solving sets of equations with matrix operations"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Dimensions of a matrix = degrees of freedom of the system \n",
+ "\n",
+ "- Given set of known outputs, $y_{1},~y_{2},~...y_{N}$ \n",
+ "- set of equations \n",
+ "- unknown inputs, $x_{1},~x_{2},~...x_{N}$, \n",
+ "- these can be written in a vector matrix format as:\n",
+ "\n",
+ "$y=Ax$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Example \n",
+ "\n",
+ "Consider a problem with 2 DOF:\n",
+ "\n",
+ "$x_{1}+3x_{2}=1$\n",
+ "\n",
+ "$2x_{1}+x_{2}=1$\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "1 & 3 \\\\\n",
+ "2 & 1 \\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",
+ "The solution for $x_{1}$ and $x_{2}$ is the intersection of two lines:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x21=[-2:2];\n",
+ "x11=1-3*x21;\n",
+ "x21=[-2:2];\n",
+ "x22=1-2*x21;\n",
+ "plot(x11,x21,x21,x22)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.40000\n",
+ " 0.20000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[1,3;2,1]; y=[1;1];\n",
+ "A\\y % matlab's Ax=y solution for x"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## 3 Degrees of Freedom\n",
+ "\n",
+ "For a $3\\times3$ matrix, the solution is the intersection of the 3 planes.\n",
+ "\n",
+ "$10x_{1}+2x_{2}+x_{3}=1$\n",
+ "\n",
+ "$2x_{1}+x_{2}+x_{3}=1$\n",
+ "\n",
+ "$x_{1}+2x_{2}+10x_{3}=1$\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "10 & 2 & 1\\\\\n",
+ "2 & 1 & 1 \\\\\n",
+ "1 & 2 & 10\\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\\\\n",
+ "x_{3} \\end{array}\\right]=\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "1 \\\\\n",
+ "1\\end{array}\\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=25;\n",
+ "x11=linspace(-2,2,N);\n",
+ "x12=linspace(-2,2,N);\n",
+ "[X11,X12]=meshgrid(x11,x12);\n",
+ "X13=1-10*X11-2*X12;\n",
+ "\n",
+ "x21=linspace(-2,2,N);\n",
+ "x22=linspace(-2,2,N);\n",
+ "[X21,X22]=meshgrid(x21,x22);\n",
+ "X23=1-2*X11-X22;\n",
+ "\n",
+ "x31=linspace(-2,2,N);\n",
+ "x32=linspace(-2,2,N);\n",
+ "[X31,X32]=meshgrid(x31,x32);\n",
+ "X33=1/10*(1-X31-2*X32);\n",
+ "\n",
+ "mesh(X11,X12,X13);\n",
+ "hold on;\n",
+ "mesh(X21,X22,X23)\n",
+ "mesh(X31,X32,X33)\n",
+ "x=[10,2, 1;2,1, 1; 1, 2, 10]\\[1;1;1];\n",
+ "plot3(x(1),x(2),x(3),'o')\n",
+ "xlabel('x1')\n",
+ "ylabel('x2')\n",
+ "zlabel('x3')\n",
+ "view(10,45)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## N>3 Degrees of Freedoms\n",
+ "\n",
+ "After 3 DOF problems, the solutions are described as *hyperplane* intersections. Which are even harder to visualize"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Gauss elimination\n",
+ "### Solving sets of equations systematically\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "2 & 1 & 1 & 1 \\\\\n",
+ "1 & 2 & 10 & 1\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Ay(2,:)-Ay(1,:)/5 = ([2 1 1 1]-1/5[10 2 1 1])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "1 & 2 & 10 & 1\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Ay(3,:)-Ay(1,:)/10 = ([1 2 10 1]-1/10[10 2 1 1])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "0 & 1.8 & 9.9 & 0.9\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Ay(3,:)-1.8\\*5/3\\*Ay(2,:) = ([0 1.8 9.9 0.9]-3\\*[0 3/5 4/5 4/5])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "0 & 0 & 7.5 & -1.5\\end{array} \n",
+ "\\right] $\n",
+ "\n",
+ "now, $7.5x_{3}=-1.5$ so $x_{3}=-\\frac{1}{5}$\n",
+ "\n",
+ "then, $3/5x_{2}+4/5(-1/5)=1$ so $x_{2}=\\frac{8}{5}$\n",
+ "\n",
+ "finally, $10x_{1}+2(8/5)+1(-\\frac{1}{5})=1$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Mass-spring example\n",
+ "\n",
+ "4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\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": 3,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 200 -100 0 0\n",
+ " -100 200 -100 0\n",
+ " 0 -100 200 -100\n",
+ " 0 0 -100 100\n",
+ "\n",
+ "y =\n",
+ "\n",
+ " 9.8100\n",
+ " 19.6200\n",
+ " 29.4300\n",
+ " 39.2400\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "k=100; % 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": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K1 =\n",
+ "\n",
+ " 200.00000 -100.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 150.00000 -100.00000 0.00000 24.52500\n",
+ " 0.00000 -100.00000 200.00000 -100.00000 29.43000\n",
+ " 0.00000 0.00000 -100.00000 100.00000 39.24000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K1=[K y];\n",
+ "K1(2,:)=K1(1,:)/2+K1(2,:)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K2 =\n",
+ "\n",
+ " 200.00000 -100.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 150.00000 -100.00000 0.00000 24.52500\n",
+ " 0.00000 0.00000 133.33333 -100.00000 45.78000\n",
+ " 0.00000 0.00000 -100.00000 100.00000 39.24000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K2=K1;\n",
+ "K2(3,:)=K1(2,:)*100/150+K1(3,:)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K2 =\n",
+ "\n",
+ " 200.00000 -100.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 150.00000 -100.00000 0.00000 24.52500\n",
+ " 0.00000 0.00000 133.33333 -100.00000 45.78000\n",
+ " 0.00000 0.00000 0.00000 25.00000 73.57500\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K2(4,:)=-K2(3,:)*K2(4,3)/K2(3,3)+K2(4,:)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x4 = 2.9430\n",
+ "x3 = 2.5506\n",
+ "x2 = 1.8639\n",
+ "x1 = 0.98100\n"
+ ]
+ }
+ ],
+ "source": [
+ "yp=K2(:,5);\n",
+ "x4=yp(4)/K2(4,4)\n",
+ "x3=(yp(3)-K2(3,4)*x4)/K2(3,3)\n",
+ "x2=(yp(2)-K2(2,3)*x3)/K2(2,2)\n",
+ "x1=(yp(1)-K2(1,2)*x2)/K2(1,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.98100\n",
+ " 1.86390\n",
+ " 2.55060\n",
+ " 2.94300\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K\\y"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Automate Gauss Elimination\n",
+ "\n",
+ "We can automate Gauss elimination with a function whose input is A and y:\n",
+ "\n",
+ "`x=GaussNaive(A,y)`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.98100\n",
+ " 1.86390\n",
+ " 2.55060\n",
+ " 2.94300\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x=GaussNaive(K,y)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Problem (Diagonal element is zero)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "If a diagonal element is 0 or very small either:\n",
+ "\n",
+ "1. no solution found\n",
+ "2. errors are introduced \n",
+ "\n",
+ "Therefore, we would want to pivot before applying Gauss elimination\n",
+ "\n",
+ "Consider:\n",
+ "\n",
+ "(a) $\\left[ \\begin{array}{cccc}\n",
+ "0 & 2 & 3 \\\\\n",
+ "4 & 6 & 7 \\\\\n",
+ "2 & -3 & 6 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\\\\n",
+ "x_{3} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "8 \\\\\n",
+ "-3 \\\\\n",
+ "5\\end{array} \\right]$\n",
+ "\n",
+ "(b) $\\left[ \\begin{array}{cccc}\n",
+ "0.0003 & 3.0000 \\\\\n",
+ "1.0000 & 1.0000 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "2.0001 \\\\\n",
+ "1.0000 \\end{array} \\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " NaN\n",
+ " NaN\n",
+ " NaN\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -5.423913\n",
+ " 0.021739\n",
+ " 2.652174\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "format short\n",
+ "Aa=[0,2,3;4,6,7;2,-3,6]; ya=[8;-3;5];\n",
+ "GaussNaive(Aa,ya)\n",
+ "Aa\\ya"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " -5.423913\n",
+ " 0.021739\n",
+ " 2.652174\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 4.00000 6.00000 7.00000 -3.00000\n",
+ " 0.00000 -6.00000 2.50000 6.50000\n",
+ " 0.00000 0.00000 3.83333 10.16667\n",
+ "\n",
+ "npivots = 2\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(Aa,ya)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.325665420556713\n",
+ " 0.666666666666667\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "format long\n",
+ "Ab=[0.3E-13,3.0000;1.0000,1.0000];yb=[2+0.1e-13;1.0000];\n",
+ "GaussNaive(Ab,yb)\n",
+ "Ab\\yb"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 1.000000000000000 1.000000000000000 1.000000000000000\n",
+ " 0.000000000000000 2.999999999999970 1.999999999999980\n",
+ "\n",
+ "npivots = 1\n",
+ "ans =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(Ab,yb)\n",
+ "Ab\\yb\n",
+ "format short"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = -3.0000\n",
+ "ans = 3.0000\n"
+ ]
+ }
+ ],
+ "source": [
+ "% determinant is (-1)^(number_of_pivots)*diagonal_elements\n",
+ "det(Ab)\n",
+ "Aug(1,1)*Aug(2,2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### Spring-Mass System again\n",
+ "Now, 4 masses are connected in series to 4 springs with $K_{1}$=10 N/m, $K_{2}$=5 N/m, \n",
+ "$K_{3}$=2 N/m \n",
+ "and $K_{4}$=1 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\n",
+ "\n",
+ "The masses have 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": 24,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 15 -5 0 0\n",
+ " -5 7 -2 0\n",
+ " 0 -2 3 -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; k2=5;k3=2;k4=1; % N/m\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": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Tridiagonal matrix\n",
+ "\n",
+ "This matrix, K, could be rewritten as 3 vectors e, f and g\n",
+ "\n",
+ "$e=\\left[ \\begin{array}{c}\n",
+ "0 \\\\\n",
+ "-5 \\\\\n",
+ "-2 \\\\\n",
+ "-1 \\end{array} \\right]$\n",
+ "\n",
+ "$f=\\left[ \\begin{array}{c}\n",
+ "15 \\\\\n",
+ "7 \\\\\n",
+ "3 \\\\\n",
+ "1 \\end{array} \\right]$\n",
+ "\n",
+ "$g=\\left[ \\begin{array}{c}\n",
+ "-5 \\\\\n",
+ "-2 \\\\\n",
+ "-1 \\\\\n",
+ "0 \\end{array} \\right]$\n",
+ "\n",
+ "Where all other components are 0 and the length of the vectors are n and the first component of e and the last component of g are zero\n",
+ "\n",
+ "`e(1)=0` \n",
+ "\n",
+ "`g(end)=0`\n",
+ "\n",
+ "No need to pivot and number of calculations reduced enormously.\n",
+ "\n",
+ "|method |Number of Floating point operations for n$\\times$n-matrix|\n",
+ "|----------------|---------|\n",
+ "| Gauss | n-cubed |\n",
+ "| Tridiagonal | n |"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 9.8100 27.4680 61.8030 101.0430\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 9.8100\n",
+ " 27.4680\n",
+ " 61.8030\n",
+ " 101.0430\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=[0;-5;-2;-1];\n",
+ "g=[-5;-2;-1;0];\n",
+ "f=[15;7;3;1];\n",
+ "Tridiag(e,f,g,y)\n",
+ "K\\y\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "% tic ... t=toc \n",
+ "% is Matlab timer used for debugging programs\n",
+ "t_GE = zeros(1,100);\n",
+ "t_GE_tridiag = zeros(1,100);\n",
+ "t_TD = zeros(1,100);\n",
+ "%for n = 1:200\n",
+ "for n=1:100\n",
+ " A = rand(n,n);\n",
+ " e = rand(n,1); e(1)=0;\n",
+ " f = rand(n,1);\n",
+ " g = rand(n,1); g(end)=0;\n",
+ " Atd=diag(f, 0) - diag(e(2:n), -1) - diag(g(1:n-1), 1);\n",
+ " b = rand(n,1);\n",
+ " tic;\n",
+ " x = GaussPivot(A,b);\n",
+ " t_GE(n) = toc;\n",
+ " tic;\n",
+ " x = A\\b;\n",
+ " t_GE_tridiag(n) = toc;\n",
+ " tic;\n",
+ " x = Tridiag(e,f,g,b);\n",
+ " t_TD(n) = toc;\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "n=1:100;\n",
+ "loglog(n,t_GE,n,t_TD,n,t_GE_tridiag)\n",
+ "legend('Gauss elim','Matlab \\','TriDiag')\n",
+ "xlabel('number of elements')\n",
+ "ylabel('time (s)')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 9.8100\n",
+ " 27.4680\n",
+ " 61.8030\n",
+ " 101.0430\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 15.00000 -5.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 5.33333 -2.00000 0.00000 22.89000\n",
+ " 0.00000 0.00000 2.25000 -1.00000 38.01375\n",
+ " 0.00000 0.00000 0.00000 0.55556 56.13500\n",
+ "\n",
+ "npivots = 0\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(K,y)"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/10_Gauss_elimination/.ipynb_checkpoints/lecture_10-checkpoint.ipynb b/10_Gauss_elimination/.ipynb_checkpoints/lecture_10-checkpoint.ipynb
new file mode 100644
index 0000000..77a6c55
--- /dev/null
+++ b/10_Gauss_elimination/.ipynb_checkpoints/lecture_10-checkpoint.ipynb
@@ -0,0 +1,1802 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Gauss Elimination\n",
+ "### Solving sets of equations with matrix operations"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Dimensions of a matrix = degrees of freedom of the system \n",
+ "\n",
+ "- Given set of known outputs, $y_{1},~y_{2},~...y_{N}$ \n",
+ "- set of equations \n",
+ "- unknown inputs, $x_{1},~x_{2},~...x_{N}$, \n",
+ "- these can be written in a vector matrix format as:\n",
+ "\n",
+ "$y=Ax$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Example \n",
+ "\n",
+ "Consider a problem with 2 DOF:\n",
+ "\n",
+ "$x_{1}+3x_{2}=1$\n",
+ "\n",
+ "$2x_{1}+x_{2}=1$\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "1 & 3 \\\\\n",
+ "2 & 1 \\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",
+ "The solution for $x_{1}$ and $x_{2}$ is the intersection of two lines:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x21=[-2:2];\n",
+ "x11=1-3*x21;\n",
+ "x21=[-2:2];\n",
+ "x22=1-2*x21;\n",
+ "plot(x11,x21,x21,x22)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.40000\n",
+ " 0.20000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[1,3;2,1]; y=[1;1];\n",
+ "A\\y % matlab's Ax=y solution for x"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## 3 Degrees of Freedom\n",
+ "\n",
+ "For a $3\\times3$ matrix, the solution is the intersection of the 3 planes.\n",
+ "\n",
+ "$10x_{1}+2x_{2}+x_{3}=1$\n",
+ "\n",
+ "$2x_{1}+x_{2}+x_{3}=1$\n",
+ "\n",
+ "$x_{1}+2x_{2}+10x_{3}=1$\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "10 & 2 & 1\\\\\n",
+ "2 & 1 & 1 \\\\\n",
+ "1 & 2 & 10\\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\\\\n",
+ "x_{3} \\end{array}\\right]=\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "1 \\\\\n",
+ "1\\end{array}\\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=25;\n",
+ "x11=linspace(-2,2,N);\n",
+ "x12=linspace(-2,2,N);\n",
+ "[X11,X12]=meshgrid(x11,x12);\n",
+ "X13=1-10*X11-2*X12;\n",
+ "\n",
+ "x21=linspace(-2,2,N);\n",
+ "x22=linspace(-2,2,N);\n",
+ "[X21,X22]=meshgrid(x21,x22);\n",
+ "X23=1-2*X11-X22;\n",
+ "\n",
+ "x31=linspace(-2,2,N);\n",
+ "x32=linspace(-2,2,N);\n",
+ "[X31,X32]=meshgrid(x31,x32);\n",
+ "X33=1/10*(1-X31-2*X32);\n",
+ "\n",
+ "mesh(X11,X12,X13);\n",
+ "hold on;\n",
+ "mesh(X21,X22,X23)\n",
+ "mesh(X31,X32,X33)\n",
+ "x=[10,2, 1;2,1, 1; 1, 2, 10]\\[1;1;1];\n",
+ "plot3(x(1),x(2),x(3),'o')\n",
+ "xlabel('x1')\n",
+ "ylabel('x2')\n",
+ "zlabel('x3')\n",
+ "view(10,45)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## N>3 Degrees of Freedoms\n",
+ "\n",
+ "After 3 DOF problems, the solutions are described as *hyperplane* intersections. Which are even harder to visualize"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Gauss elimination\n",
+ "### Solving sets of equations systematically\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "2 & 1 & 1 & 1 \\\\\n",
+ "1 & 2 & 10 & 1\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Ay(2,:)-Ay(1,:)/5 = ([2 1 1 1]-1/5[10 2 1 1])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "1 & 2 & 10 & 1\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Ay(3,:)-Ay(1,:)/10 = ([1 2 10 1]-1/10[10 2 1 1])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "0 & 1.8 & 9.9 & 0.9\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Ay(3,:)-1.8\\*5/3\\*Ay(2,:) = ([0 1.8 9.9 0.9]-3\\*[0 3/5 4/5 4/5])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "0 & 0 & 7.5 & -1.5\\end{array} \n",
+ "\\right] $\n",
+ "\n",
+ "now, $7.5x_{3}=-1.5$ so $x_{3}=-\\frac{1}{5}$\n",
+ "\n",
+ "then, $3/5x_{2}+4/5(-1/5)=1$ so $x_{2}=\\frac{8}{5}$\n",
+ "\n",
+ "finally, $10x_{1}+2(8/5)+1(-\\frac{1}{5})=1$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Mass-spring example\n",
+ "\n",
+ "4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K1 =\n",
+ "\n",
+ " 20.00000 -10.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 15.00000 -10.00000 0.00000 24.52500\n",
+ " 0.00000 -10.00000 20.00000 -10.00000 29.43000\n",
+ " 0.00000 0.00000 -10.00000 10.00000 39.24000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K1=[K y];\n",
+ "K1(2,:)=K1(1,:)/2+K1(2,:)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K2 =\n",
+ "\n",
+ " 20.00000 -10.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 15.00000 -10.00000 0.00000 24.52500\n",
+ " 0.00000 0.00000 13.33333 -10.00000 45.78000\n",
+ " 0.00000 0.00000 -10.00000 10.00000 39.24000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K2=K1;\n",
+ "K2(3,:)=K1(2,:)*2/3+K1(3,:)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K2 =\n",
+ "\n",
+ " 20.00000 -10.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 15.00000 -10.00000 0.00000 24.52500\n",
+ " 0.00000 0.00000 13.33333 -10.00000 45.78000\n",
+ " 0.00000 0.00000 0.00000 2.50000 73.57500\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K2(4,:)=-K2(3,:)*K2(4,3)/K2(3,3)+K2(4,:)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x4 = 29.430\n",
+ "x3 = 25.506\n",
+ "x2 = 18.639\n",
+ "x1 = 9.8100\n"
+ ]
+ }
+ ],
+ "source": [
+ "yp=K2(:,5);\n",
+ "x4=yp(4)/K2(4,4)\n",
+ "x3=(yp(3)+10*x4)/K2(3,3)\n",
+ "x2=(yp(2)+10*x3)/K2(2,2)\n",
+ "x1=(yp(1)+10*x2)/K2(1,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 9.8100\n",
+ " 18.6390\n",
+ " 25.5060\n",
+ " 29.4300\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K\\y"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Automate Gauss Elimination\n",
+ "\n",
+ "We can automate Gauss elimination with a function whose input is A and y:\n",
+ "\n",
+ "`x=GaussNaive(A,y)`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 9.8100\n",
+ " 18.6390\n",
+ " 25.5060\n",
+ " 29.4300\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x=GaussNaive(K,y)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Problem (Diagonal element is zero)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "If a diagonal element is 0 or very small either:\n",
+ "\n",
+ "1. no solution found\n",
+ "2. errors are introduced \n",
+ "\n",
+ "Therefore, we would want to pivot before applying Gauss elimination\n",
+ "\n",
+ "Consider:\n",
+ "\n",
+ "(a) $\\left[ \\begin{array}{cccc}\n",
+ "0 & 2 & 3 \\\\\n",
+ "4 & 6 & 7 \\\\\n",
+ "2 & -3 & 6 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\\\\n",
+ "x_{3} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "8 \\\\\n",
+ "-3 \\\\\n",
+ "5\\end{array} \\right]$\n",
+ "\n",
+ "(b) $\\left[ \\begin{array}{cccc}\n",
+ "0.0003 & 3.0000 \\\\\n",
+ "1.0000 & 1.0000 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "2.0001 \\\\\n",
+ "1.0000 \\end{array} \\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " NaN\n",
+ " NaN\n",
+ " NaN\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -5.423913\n",
+ " 0.021739\n",
+ " 2.652174\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "format short\n",
+ "Aa=[0,2,3;4,6,7;2,-3,6]; ya=[8;-3;5];\n",
+ "GaussNaive(Aa,ya)\n",
+ "Aa\\ya"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " -5.423913\n",
+ " 0.021739\n",
+ " 2.652174\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 4.00000 6.00000 7.00000 -3.00000\n",
+ " 0.00000 -6.00000 2.50000 6.50000\n",
+ " 0.00000 0.00000 3.83333 10.16667\n",
+ "\n",
+ "npivots = 2\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(Aa,ya)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.325665420556713\n",
+ " 0.666666666666667\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "format long\n",
+ "Ab=[0.3E-13,3.0000;1.0000,1.0000];yb=[2+0.1e-13;1.0000];\n",
+ "GaussNaive(Ab,yb)\n",
+ "Ab\\yb"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 1.000000000000000 1.000000000000000 1.000000000000000\n",
+ " 0.000000000000000 2.999999999999970 1.999999999999980\n",
+ "\n",
+ "npivots = 1\n",
+ "ans =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(Ab,yb)\n",
+ "Ab\\yb\n",
+ "format short"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = -3.0000\n",
+ "ans = 3.0000\n"
+ ]
+ }
+ ],
+ "source": [
+ "% determinant is (-1)^(number_of_pivots)*diagonal_elements\n",
+ "det(Ab)\n",
+ "Aug(1,1)*Aug(2,2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### Spring-Mass System again\n",
+ "Now, 4 masses are connected in series to 4 springs with $K_{1}$=10 N/m, $K_{2}$=5 N/m, \n",
+ "$K_{3}$=2 N/m \n",
+ "and $K_{4}$=1 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\n",
+ "\n",
+ "The masses have 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": 24,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 15 -5 0 0\n",
+ " -5 7 -2 0\n",
+ " 0 -2 3 -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; k2=5;k3=2;k4=1; % N/m\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": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Tridiagonal matrix\n",
+ "\n",
+ "This matrix, K, could be rewritten as 3 vectors e, f and g\n",
+ "\n",
+ "$e=\\left[ \\begin{array}{c}\n",
+ "0 \\\\\n",
+ "-5 \\\\\n",
+ "-2 \\\\\n",
+ "-1 \\end{array} \\right]$\n",
+ "\n",
+ "$f=\\left[ \\begin{array}{c}\n",
+ "15 \\\\\n",
+ "7 \\\\\n",
+ "3 \\\\\n",
+ "1 \\end{array} \\right]$\n",
+ "\n",
+ "$g=\\left[ \\begin{array}{c}\n",
+ "-5 \\\\\n",
+ "-2 \\\\\n",
+ "-1 \\\\\n",
+ "0 \\end{array} \\right]$\n",
+ "\n",
+ "Where all other components are 0 and the length of the vectors are n and the first component of e and the last component of g are zero\n",
+ "\n",
+ "`e(1)=0` \n",
+ "\n",
+ "`g(end)=0`\n",
+ "\n",
+ "No need to pivot and number of calculations reduced enormously.\n",
+ "\n",
+ "|method |Number of Floating point operations for n$\\times$n-matrix|\n",
+ "|----------------|---------|\n",
+ "| Gauss | n-cubed |\n",
+ "| Tridiagonal | n |"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 9.8100 27.4680 61.8030 101.0430\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 9.8100\n",
+ " 27.4680\n",
+ " 61.8030\n",
+ " 101.0430\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=[0;-5;-2;-1];\n",
+ "g=[-5;-2;-1;0];\n",
+ "f=[15;7;3;1];\n",
+ "Tridiag(e,f,g,y)\n",
+ "K\\y\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "% tic ... t=toc \n",
+ "% is Matlab timer used for debugging programs\n",
+ "t_GE = zeros(1,100);\n",
+ "t_GE_tridiag = zeros(1,100);\n",
+ "t_TD = zeros(1,100);\n",
+ "%for n = 1:200\n",
+ "for n=1:100\n",
+ " A = rand(n,n);\n",
+ " e = rand(n,1); e(1)=0;\n",
+ " f = rand(n,1);\n",
+ " g = rand(n,1); g(end)=0;\n",
+ " Atd=diag(f, 0) - diag(e(2:n), -1) - diag(g(1:n-1), 1);\n",
+ " b = rand(n,1);\n",
+ " tic;\n",
+ " x = GaussPivot(A,b);\n",
+ " t_GE(n) = toc;\n",
+ " tic;\n",
+ " x = A\\b;\n",
+ " t_GE_tridiag(n) = toc;\n",
+ " tic;\n",
+ " x = Tridiag(e,f,g,b);\n",
+ " t_TD(n) = toc;\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "n=1:100;\n",
+ "loglog(n,t_GE,n,t_TD,n,t_GE_tridiag)\n",
+ "legend('Gauss elim','Matlab \\','TriDiag')\n",
+ "xlabel('number of elements')\n",
+ "ylabel('time (s)')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 9.8100\n",
+ " 27.4680\n",
+ " 61.8030\n",
+ " 101.0430\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 15.00000 -5.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 5.33333 -2.00000 0.00000 22.89000\n",
+ " 0.00000 0.00000 2.25000 -1.00000 38.01375\n",
+ " 0.00000 0.00000 0.00000 0.55556 56.13500\n",
+ "\n",
+ "npivots = 0\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(K,y)"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/10_Gauss_elimination/10_Gauss_elimination.ipynb b/10_Gauss_elimination/10_Gauss_elimination.ipynb
new file mode 100644
index 0000000..780de50
--- /dev/null
+++ b/10_Gauss_elimination/10_Gauss_elimination.ipynb
@@ -0,0 +1,1786 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Gauss Elimination\n",
+ "### Solving sets of equations with matrix operations"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Dimensions of unknown vector = degrees of freedom of the system \n",
+ "\n",
+ "- Given set of known outputs, $y_{1},~y_{2},~...y_{N}$ \n",
+ "- set of equations \n",
+ "- unknown inputs, $x_{1},~x_{2},~...x_{N}$, \n",
+ "- these can be written in a vector matrix format as:\n",
+ "\n",
+ "$y=Ax$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Example \n",
+ "\n",
+ "Consider a problem with 2 DOF:\n",
+ "\n",
+ "$x_{1}+3x_{2}=1$\n",
+ "\n",
+ "$2x_{1}+x_{2}=1$\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "1 & 3 \\\\\n",
+ "2 & 1 \\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",
+ "The solution for $x_{1}$ and $x_{2}$ is the intersection of two lines:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x21=[-2:2];\n",
+ "x11=1-3*x21;\n",
+ "x21=[-2:2];\n",
+ "x22=1-2*x21;\n",
+ "plot(x11,x21,x21,x22)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.40000\n",
+ " 0.20000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[1,3;2,1]; y=[1;1];\n",
+ "A\\y % matlab's Ax=y solution for x"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## 3 Degrees of Freedom\n",
+ "\n",
+ "For a $3\\times3$ matrix, the solution is the intersection of the 3 planes.\n",
+ "\n",
+ "$10x_{1}+2x_{2}+x_{3}=1$\n",
+ "\n",
+ "$2x_{1}+x_{2}+x_{3}=1$\n",
+ "\n",
+ "$x_{1}+2x_{2}+10x_{3}=1$\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "10 & 2 & 1\\\\\n",
+ "2 & 1 & 1 \\\\\n",
+ "1 & 2 & 10\\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\\\\n",
+ "x_{3} \\end{array}\\right]=\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "1 \\\\\n",
+ "1\\end{array}\\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "mesh(X11,X12,X13);\n",
+ "hold on;\n",
+ "mesh(X21,X22,X23)\n",
+ "mesh(X31,X32,X33)\n",
+ "x=[10,2, 1;2,1, 1; 1, 2, 10]\\[1;1;1];\n",
+ "plot3(x(1),x(2),x(3),'o')\n",
+ "xlabel('x1')\n",
+ "ylabel('x2')\n",
+ "zlabel('x3')\n",
+ "view(10,45)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## N>3 Degrees of Freedoms\n",
+ "\n",
+ "After 3 DOF problems, the solutions are described as *hyperplane* intersections. Which are even harder to visualize"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Gauss elimination\n",
+ "### Solving sets of equations systematically\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "2 & 1 & 1 & 1 \\\\\n",
+ "1 & 2 & 10 & 1\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Ay(2,:)-Ay(1,:)/5 = ([2 1 1 1]-1/5[10 2 1 1])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "1 & 2 & 10 & 1\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Ay(3,:)-Ay(1,:)/10 = ([1 2 10 1]-1/10[10 2 1 1])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "0 & 1.8 & 9.9 & 0.9\\end{array} \n",
+ "\\right] $"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Ay(3,:)-1.8\\*5/3\\*Ay(2,:) = ([0 1.8 9.9 0.9]-3\\*[0 3/5 4/5 4/5])\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc|c}\n",
+ " & A & & y \\\\\n",
+ "10 & 2 & 1 & 1\\\\\n",
+ "0 & 3/5 & 4/5 & 4/5 \\\\\n",
+ "0 & 0 & 7.5 & -1.5\\end{array} \n",
+ "\\right] $\n",
+ "\n",
+ "now, $7.5x_{3}=-1.5$ so $x_{3}=-\\frac{1}{5}$\n",
+ "\n",
+ "then, $3/5x_{2}+4/5(-1/5)=1$ so $x_{2}=\\frac{8}{5}$\n",
+ "\n",
+ "finally, $10x_{1}+2(8/5)+1(-\\frac{1}{5})=1$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Mass-spring example\n",
+ "\n",
+ "4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\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": 3,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 200 -100 0 0\n",
+ " -100 200 -100 0\n",
+ " 0 -100 200 -100\n",
+ " 0 0 -100 100\n",
+ "\n",
+ "y =\n",
+ "\n",
+ " 9.8100\n",
+ " 19.6200\n",
+ " 29.4300\n",
+ " 39.2400\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "k=100; % 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": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K1 =\n",
+ "\n",
+ " 200.00000 -100.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 150.00000 -100.00000 0.00000 24.52500\n",
+ " 0.00000 -100.00000 200.00000 -100.00000 29.43000\n",
+ " 0.00000 0.00000 -100.00000 100.00000 39.24000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K1=[K y];\n",
+ "K1(2,:)=K1(1,:)/2+K1(2,:)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K2 =\n",
+ "\n",
+ " 200.00000 -100.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 150.00000 -100.00000 0.00000 24.52500\n",
+ " 0.00000 0.00000 133.33333 -100.00000 45.78000\n",
+ " 0.00000 0.00000 -100.00000 100.00000 39.24000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K2=K1;\n",
+ "K2(3,:)=K1(2,:)*100/150+K1(3,:)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K2 =\n",
+ "\n",
+ " 200.00000 -100.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 150.00000 -100.00000 0.00000 24.52500\n",
+ " 0.00000 0.00000 133.33333 -100.00000 45.78000\n",
+ " 0.00000 0.00000 0.00000 25.00000 73.57500\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K2(4,:)=-K2(3,:)*K2(4,3)/K2(3,3)+K2(4,:)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x4 = 2.9430\n",
+ "x3 = 2.5506\n",
+ "x2 = 1.8639\n",
+ "x1 = 0.98100\n"
+ ]
+ }
+ ],
+ "source": [
+ "yp=K2(:,5);\n",
+ "x4=yp(4)/K2(4,4)\n",
+ "x3=(yp(3)-K2(3,4)*x4)/K2(3,3)\n",
+ "x2=(yp(2)-K2(2,3)*x3)/K2(2,2)\n",
+ "x1=(yp(1)-K2(1,2)*x2)/K2(1,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.98100\n",
+ " 1.86390\n",
+ " 2.55060\n",
+ " 2.94300\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K\\y"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Automate Gauss Elimination\n",
+ "\n",
+ "We can automate Gauss elimination with a function whose input is A and y:\n",
+ "\n",
+ "`x=GaussNaive(A,y)`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.98100\n",
+ " 1.86390\n",
+ " 2.55060\n",
+ " 2.94300\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x=GaussNaive(K,y)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Problem (Diagonal element is zero)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "If a diagonal element is 0 or very small either:\n",
+ "\n",
+ "1. no solution found\n",
+ "2. errors are introduced \n",
+ "\n",
+ "Therefore, we would want to pivot before applying Gauss elimination\n",
+ "\n",
+ "Consider:\n",
+ "\n",
+ "(a) $\\left[ \\begin{array}{cccc}\n",
+ "0 & 2 & 3 \\\\\n",
+ "4 & 6 & 7 \\\\\n",
+ "2 & -3 & 6 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\\\\n",
+ "x_{3} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "8 \\\\\n",
+ "-3 \\\\\n",
+ "5\\end{array} \\right]$\n",
+ "\n",
+ "(b) $\\left[ \\begin{array}{cccc}\n",
+ "0.0003 & 3.0000 \\\\\n",
+ "1.0000 & 1.0000 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "2.0001 \\\\\n",
+ "1.0000 \\end{array} \\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " NaN\n",
+ " NaN\n",
+ " NaN\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -5.423913\n",
+ " 0.021739\n",
+ " 2.652174\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "format short\n",
+ "Aa=[0,2,3;4,6,7;2,-3,6]; ya=[8;-3;5];\n",
+ "GaussNaive(Aa,ya)\n",
+ "Aa\\ya"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " -5.423913\n",
+ " 0.021739\n",
+ " 2.652174\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 4.00000 6.00000 7.00000 -3.00000\n",
+ " 0.00000 -6.00000 2.50000 6.50000\n",
+ " 0.00000 0.00000 3.83333 10.16667\n",
+ "\n",
+ "npivots = 2\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(Aa,ya)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 0.325665420556713\n",
+ " 0.666666666666667\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "format long\n",
+ "Ab=[0.3E-13,3.0000;1.0000,1.0000];yb=[2+0.1e-13;1.0000];\n",
+ "GaussNaive(Ab,yb)\n",
+ "Ab\\yb"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 1.000000000000000 1.000000000000000 1.000000000000000\n",
+ " 0.000000000000000 2.999999999999970 1.999999999999980\n",
+ "\n",
+ "npivots = 1\n",
+ "ans =\n",
+ "\n",
+ " 0.333333333333333\n",
+ " 0.666666666666667\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(Ab,yb)\n",
+ "Ab\\yb\n",
+ "format short"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = -3.0000\n",
+ "ans = 3.0000\n"
+ ]
+ }
+ ],
+ "source": [
+ "% determinant is (-1)^(number_of_pivots)*diagonal_elements\n",
+ "det(Ab)\n",
+ "Aug(1,1)*Aug(2,2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### Spring-Mass System again\n",
+ "Now, 4 masses are connected in series to 4 springs with $K_{1}$=10 N/m, $K_{2}$=5 N/m, \n",
+ "$K_{3}$=2 N/m \n",
+ "and $K_{4}$=1 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\n",
+ "\n",
+ "The masses have 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": 24,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 15 -5 0 0\n",
+ " -5 7 -2 0\n",
+ " 0 -2 3 -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; k2=5;k3=2;k4=1; % N/m\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": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Tridiagonal matrix\n",
+ "\n",
+ "This matrix, K, could be rewritten as 3 vectors e, f and g\n",
+ "\n",
+ "$e=\\left[ \\begin{array}{c}\n",
+ "0 \\\\\n",
+ "-5 \\\\\n",
+ "-2 \\\\\n",
+ "-1 \\end{array} \\right]$\n",
+ "\n",
+ "$f=\\left[ \\begin{array}{c}\n",
+ "15 \\\\\n",
+ "7 \\\\\n",
+ "3 \\\\\n",
+ "1 \\end{array} \\right]$\n",
+ "\n",
+ "$g=\\left[ \\begin{array}{c}\n",
+ "-5 \\\\\n",
+ "-2 \\\\\n",
+ "-1 \\\\\n",
+ "0 \\end{array} \\right]$\n",
+ "\n",
+ "Where all other components are 0 and the length of the vectors are n and the first component of e and the last component of g are zero\n",
+ "\n",
+ "`e(1)=0` \n",
+ "\n",
+ "`g(end)=0`\n",
+ "\n",
+ "No need to pivot and number of calculations reduced enormously.\n",
+ "\n",
+ "|method |Number of Floating point operations for n$\\times$n-matrix|\n",
+ "|----------------|---------|\n",
+ "| Gauss | n-cubed |\n",
+ "| Tridiagonal | n |"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 9.8100 27.4680 61.8030 101.0430\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 9.8100\n",
+ " 27.4680\n",
+ " 61.8030\n",
+ " 101.0430\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=[0;-5;-2;-1];\n",
+ "g=[-5;-2;-1;0];\n",
+ "f=[15;7;3;1];\n",
+ "Tridiag(e,f,g,y)\n",
+ "K\\y\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "% tic ... t=toc \n",
+ "% is Matlab timer used for debugging programs\n",
+ "t_GE = zeros(1,100);\n",
+ "t_GE_tridiag = zeros(1,100);\n",
+ "t_TD = zeros(1,100);\n",
+ "%for n = 1:200\n",
+ "for n=1:100\n",
+ " A = rand(n,n);\n",
+ " e = rand(n,1); e(1)=0;\n",
+ " f = rand(n,1);\n",
+ " g = rand(n,1); g(end)=0;\n",
+ " Atd=diag(f, 0) - diag(e(2:n), -1) - diag(g(1:n-1), 1);\n",
+ " b = rand(n,1);\n",
+ " tic;\n",
+ " x = GaussPivot(A,b);\n",
+ " t_GE(n) = toc;\n",
+ " tic;\n",
+ " x = A\\b;\n",
+ " t_GE_tridiag(n) = toc;\n",
+ " tic;\n",
+ " x = Tridiag(e,f,g,b);\n",
+ " t_TD(n) = toc;\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "n=1:100;\n",
+ "loglog(n,t_GE,n,t_TD,n,t_GE_tridiag)\n",
+ "legend('Gauss elim','Matlab \\','TriDiag')\n",
+ "xlabel('number of elements')\n",
+ "ylabel('time (s)')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 9.8100\n",
+ " 27.4680\n",
+ " 61.8030\n",
+ " 101.0430\n",
+ "\n",
+ "Aug =\n",
+ "\n",
+ " 15.00000 -5.00000 0.00000 0.00000 9.81000\n",
+ " 0.00000 5.33333 -2.00000 0.00000 22.89000\n",
+ " 0.00000 0.00000 2.25000 -1.00000 38.01375\n",
+ " 0.00000 0.00000 0.00000 0.55556 56.13500\n",
+ "\n",
+ "npivots = 0\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,Aug,npivots]=GaussPivot(K,y)"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/10_Gauss_elimination/GaussNaive.m b/10_Gauss_elimination/GaussNaive.m
new file mode 100644
index 0000000..d8c60d3
--- /dev/null
+++ b/10_Gauss_elimination/GaussNaive.m
@@ -0,0 +1,25 @@
+function [x,Aug] = GaussNaive(A,y)
+% 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+1;
+Aug = [A y];
+% forward elimination
+for k = 1:n-1
+ for i = k+1:n
+ factor = Aug(i,k)/Aug(k,k);
+ Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(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/10_Gauss_elimination/GaussPivot.m b/10_Gauss_elimination/GaussPivot.m
new file mode 100644
index 0000000..20df9e8
--- /dev/null
+++ b/10_Gauss_elimination/GaussPivot.m
@@ -0,0 +1,33 @@
+function [x,Aug,npivots] = GaussPivot(A,b)
+% GaussPivot: Gauss elimination pivoting
+% x = GaussPivot(A,b): Gauss elimination with pivoting.
+% input:
+% A = coefficient matrix
+% b = right hand side vector
+% output:
+% x = solution vector
+[m,n]=size(A);
+if m~=n, error('Matrix A must be square'); end
+nb=n+1;
+Aug=[A b];
+npivots=0; % initially no pivots used
+% forward elimination
+for k = 1:n-1
+ % partial pivoting
+ [big,i]=max(abs(Aug(k:n,k)));
+ ipr=i+k-1;
+ if ipr~=k
+ npivots=npivots+1; % if the max is not the current index ipr, pivot count
+ Aug([k,ipr],:)=Aug([ipr,k],:);
+ end
+ for i = k+1:n
+ factor=Aug(i,k)/Aug(k,k);
+ Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(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/10_Gauss_elimination/Tridiag.m b/10_Gauss_elimination/Tridiag.m
new file mode 100644
index 0000000..ac4ec9b
--- /dev/null
+++ b/10_Gauss_elimination/Tridiag.m
@@ -0,0 +1,22 @@
+function x = Tridiag(e,f,g,r)
+% Tridiag: Tridiagonal equation solver banded system
+% x = Tridiag(e,f,g,r): Tridiagonal system solver.
+% input:
+% e = subdiagonal vector
+% f = diagonal vector
+% g = superdiagonal vector
+% r = right hand side vector
+% output:
+% x = solution vector
+n=length(f);
+% forward elimination
+for k = 2:n
+ factor = e(k)/f(k-1);
+ f(k) = f(k) - factor*g(k-1);
+ r(k) = r(k) - factor*r(k-1);
+end
+% back substitution
+x(n) = r(n)/f(n);
+for k = n-1:-1:1
+ x(k) = (r(k)-g(k)*x(k+1))/f(k);
+end
diff --git a/10_Gauss_elimination/nohup.out b/10_Gauss_elimination/nohup.out
new file mode 100644
index 0000000..e69de29
diff --git a/10_Gauss_elimination/octave-workspace b/10_Gauss_elimination/octave-workspace
new file mode 100644
index 0000000..08dd0b6
Binary files /dev/null and b/10_Gauss_elimination/octave-workspace differ
diff --git "a/10_Gauss_elimination/\340\004'\001" "b/10_Gauss_elimination/\340\004'\001"
new file mode 100644
index 0000000..e69de29
diff --git "a/10_Gauss_elimination/\360\206P\001" "b/10_Gauss_elimination/\360\206P\001"
new file mode 100644
index 0000000..e69de29
diff --git a/11_LU-decomposition/.ipynb_checkpoints/11_LU-decomposition-checkpoint.ipynb b/11_LU-decomposition/.ipynb_checkpoints/11_LU-decomposition-checkpoint.ipynb
new file mode 100644
index 0000000..7f37020
--- /dev/null
+++ b/11_LU-decomposition/.ipynb_checkpoints/11_LU-decomposition-checkpoint.ipynb
@@ -0,0 +1,760 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# LU Decomposition\n",
+ "### efficient storage of matrices for solutions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Consider equation set:\n",
+ "\n",
+ "$y=Ax$\n",
+ "\n",
+ "Assume that we can perform Gauss elimination and achieve this formula:\n",
+ "\n",
+ "$Ux=d$ \n",
+ "\n",
+ "$U$: upper triangular matrix derived from Gauss elimination \n",
+ "\n",
+ "$d$: is the set of dependent variables after Gauss elimination. \n",
+ "\n",
+ "Assume lower triangular matrix, $L$, exists \n",
+ "\n",
+ "- ones on the diagonal \n",
+ "- same dimensions as $U$ \n",
+ "- following is true:\n",
+ "\n",
+ "$L(Ux-d)=Ax-y=0$\n",
+ "\n",
+ "Now, $Ax=LUx$, so $A=LU$, and $y=Ld$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "$2x_{1}+x_{2}=1$\n",
+ "\n",
+ "$x_{1}+3x_{2}=1$\n",
+ "\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "2 & 1 \\\\\n",
+ "1 & 3 \\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\end{array}\\right]=\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "1\\end{array}\\right]$\n",
+ "\n",
+ "f21=0.5\n",
+ "\n",
+ "A(2,1)=1-1 = 0 \n",
+ "\n",
+ "A(2,2)=3-0.5=2.5\n",
+ "\n",
+ "y(2)=1-0.5=0.5\n",
+ "\n",
+ "$L(Ux-d)=\n",
+ "\\left[ \\begin{array}{cc}\n",
+ "1 & 0 \\\\\n",
+ "0.5 & 1 \\end{array} \\right]\n",
+ "\\left(\\left[ \\begin{array}{cc}\n",
+ "2 & 1 \\\\\n",
+ "0 & 2.5 \\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\end{array}\\right]-\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "0.5\\end{array}\\right]\\right)=0$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n",
+ "\n",
+ " -- Built-in Function: [L, U] = lu (A)\n",
+ " -- Built-in Function: [L, U, P] = lu (A)\n",
+ " -- Built-in Function: [L, U, P, Q] = lu (S)\n",
+ " -- Built-in Function: [L, U, P, Q, R] = lu (S)\n",
+ " -- Built-in Function: [...] = lu (S, THRES)\n",
+ " -- Built-in Function: Y = lu (...)\n",
+ " -- Built-in Function: [...] = lu (..., \"vector\")\n",
+ " Compute the LU decomposition of A.\n",
+ "\n",
+ " If A is full subroutines from LAPACK are used and if A is sparse\n",
+ " then UMFPACK is used.\n",
+ "\n",
+ " The result is returned in a permuted form, according to the\n",
+ " optional return value P. For example, given the matrix 'a = [1, 2;\n",
+ " 3, 4]',\n",
+ "\n",
+ " [l, u, p] = lu (A)\n",
+ "\n",
+ " returns\n",
+ "\n",
+ " l =\n",
+ "\n",
+ " 1.00000 0.00000\n",
+ " 0.33333 1.00000\n",
+ "\n",
+ " u =\n",
+ "\n",
+ " 3.00000 4.00000\n",
+ " 0.00000 0.66667\n",
+ "\n",
+ " p =\n",
+ "\n",
+ " 0 1\n",
+ " 1 0\n",
+ "\n",
+ " The matrix is not required to be square.\n",
+ "\n",
+ " When called with two or three output arguments and a spare input\n",
+ " matrix, 'lu' does not attempt to perform sparsity preserving column\n",
+ " permutations. Called with a fourth output argument, the sparsity\n",
+ " preserving column transformation Q is returned, such that 'P * A *\n",
+ " Q = L * U'.\n",
+ "\n",
+ " Called with a fifth output argument and a sparse input matrix, 'lu'\n",
+ " attempts to use a scaling factor R on the input matrix such that 'P\n",
+ " * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n",
+ " stable factorization.\n",
+ "\n",
+ " An additional input argument THRES, that defines the pivoting\n",
+ " threshold can be given. THRES can be a scalar, in which case it\n",
+ " defines the UMFPACK pivoting tolerance for both symmetric and\n",
+ " unsymmetric cases. If THRES is a 2-element vector, then the first\n",
+ " element defines the pivoting tolerance for the unsymmetric UMFPACK\n",
+ " pivoting strategy and the second for the symmetric strategy. By\n",
+ " default, the values defined by 'spparms' are used ([0.1, 0.001]).\n",
+ "\n",
+ " Given the string argument \"vector\", 'lu' returns the values of P\n",
+ " and Q as vector values, such that for full matrix, 'A (P,:) = L *\n",
+ " U', and 'R(P,:) * A (:, Q) = L * U'.\n",
+ "\n",
+ " With two output arguments, returns the permuted forms of the upper\n",
+ " and lower triangular matrices, such that 'A = L * U'. With one\n",
+ " output argument Y, then the matrix returned by the LAPACK routines\n",
+ " is returned. If the input matrix is sparse then the matrix L is\n",
+ " embedded into U to give a return value similar to the full case.\n",
+ " For both full and sparse matrices, 'lu' loses the permutation\n",
+ " information.\n",
+ "\n",
+ " See also: luupdate, ilu, chol, hess, qr, qz, schur, svd.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help lu"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% time LU solution vs backslash\n",
+ "t_lu=zeros(100,1);\n",
+ "t_bs=zeros(100,1);\n",
+ "for N=1:100\n",
+ " A=rand(N,N);\n",
+ " y=rand(N,1);\n",
+ " [L,U,P]=lu(A);\n",
+ "\n",
+ " tic; d=inv(L)*y; x=inv(U)*d; t_lu(N)=toc;\n",
+ "\n",
+ " tic; x=inv(A)*y; t_bs(N)=toc;\n",
+ "end\n",
+ "plot([1:100],t_lu,[1:100],t_bs) \n",
+ "legend('LU decomp','Octave \\\\')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Spring-mass system\n",
+ "\n",
+ "4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "This matrix, K, is symmetric. \n",
+ "\n",
+ "`K(i,j)==K(j,i)`\n",
+ "\n",
+ "Now we can use,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "## 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": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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))"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/11_LU-decomposition/.ipynb_checkpoints/lecture_11-checkpoint.ipynb b/11_LU-decomposition/.ipynb_checkpoints/lecture_11-checkpoint.ipynb
new file mode 100644
index 0000000..7f37020
--- /dev/null
+++ b/11_LU-decomposition/.ipynb_checkpoints/lecture_11-checkpoint.ipynb
@@ -0,0 +1,760 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# LU Decomposition\n",
+ "### efficient storage of matrices for solutions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Consider equation set:\n",
+ "\n",
+ "$y=Ax$\n",
+ "\n",
+ "Assume that we can perform Gauss elimination and achieve this formula:\n",
+ "\n",
+ "$Ux=d$ \n",
+ "\n",
+ "$U$: upper triangular matrix derived from Gauss elimination \n",
+ "\n",
+ "$d$: is the set of dependent variables after Gauss elimination. \n",
+ "\n",
+ "Assume lower triangular matrix, $L$, exists \n",
+ "\n",
+ "- ones on the diagonal \n",
+ "- same dimensions as $U$ \n",
+ "- following is true:\n",
+ "\n",
+ "$L(Ux-d)=Ax-y=0$\n",
+ "\n",
+ "Now, $Ax=LUx$, so $A=LU$, and $y=Ld$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "$2x_{1}+x_{2}=1$\n",
+ "\n",
+ "$x_{1}+3x_{2}=1$\n",
+ "\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "2 & 1 \\\\\n",
+ "1 & 3 \\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\end{array}\\right]=\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "1\\end{array}\\right]$\n",
+ "\n",
+ "f21=0.5\n",
+ "\n",
+ "A(2,1)=1-1 = 0 \n",
+ "\n",
+ "A(2,2)=3-0.5=2.5\n",
+ "\n",
+ "y(2)=1-0.5=0.5\n",
+ "\n",
+ "$L(Ux-d)=\n",
+ "\\left[ \\begin{array}{cc}\n",
+ "1 & 0 \\\\\n",
+ "0.5 & 1 \\end{array} \\right]\n",
+ "\\left(\\left[ \\begin{array}{cc}\n",
+ "2 & 1 \\\\\n",
+ "0 & 2.5 \\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\end{array}\\right]-\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "0.5\\end{array}\\right]\\right)=0$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n",
+ "\n",
+ " -- Built-in Function: [L, U] = lu (A)\n",
+ " -- Built-in Function: [L, U, P] = lu (A)\n",
+ " -- Built-in Function: [L, U, P, Q] = lu (S)\n",
+ " -- Built-in Function: [L, U, P, Q, R] = lu (S)\n",
+ " -- Built-in Function: [...] = lu (S, THRES)\n",
+ " -- Built-in Function: Y = lu (...)\n",
+ " -- Built-in Function: [...] = lu (..., \"vector\")\n",
+ " Compute the LU decomposition of A.\n",
+ "\n",
+ " If A is full subroutines from LAPACK are used and if A is sparse\n",
+ " then UMFPACK is used.\n",
+ "\n",
+ " The result is returned in a permuted form, according to the\n",
+ " optional return value P. For example, given the matrix 'a = [1, 2;\n",
+ " 3, 4]',\n",
+ "\n",
+ " [l, u, p] = lu (A)\n",
+ "\n",
+ " returns\n",
+ "\n",
+ " l =\n",
+ "\n",
+ " 1.00000 0.00000\n",
+ " 0.33333 1.00000\n",
+ "\n",
+ " u =\n",
+ "\n",
+ " 3.00000 4.00000\n",
+ " 0.00000 0.66667\n",
+ "\n",
+ " p =\n",
+ "\n",
+ " 0 1\n",
+ " 1 0\n",
+ "\n",
+ " The matrix is not required to be square.\n",
+ "\n",
+ " When called with two or three output arguments and a spare input\n",
+ " matrix, 'lu' does not attempt to perform sparsity preserving column\n",
+ " permutations. Called with a fourth output argument, the sparsity\n",
+ " preserving column transformation Q is returned, such that 'P * A *\n",
+ " Q = L * U'.\n",
+ "\n",
+ " Called with a fifth output argument and a sparse input matrix, 'lu'\n",
+ " attempts to use a scaling factor R on the input matrix such that 'P\n",
+ " * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n",
+ " stable factorization.\n",
+ "\n",
+ " An additional input argument THRES, that defines the pivoting\n",
+ " threshold can be given. THRES can be a scalar, in which case it\n",
+ " defines the UMFPACK pivoting tolerance for both symmetric and\n",
+ " unsymmetric cases. If THRES is a 2-element vector, then the first\n",
+ " element defines the pivoting tolerance for the unsymmetric UMFPACK\n",
+ " pivoting strategy and the second for the symmetric strategy. By\n",
+ " default, the values defined by 'spparms' are used ([0.1, 0.001]).\n",
+ "\n",
+ " Given the string argument \"vector\", 'lu' returns the values of P\n",
+ " and Q as vector values, such that for full matrix, 'A (P,:) = L *\n",
+ " U', and 'R(P,:) * A (:, Q) = L * U'.\n",
+ "\n",
+ " With two output arguments, returns the permuted forms of the upper\n",
+ " and lower triangular matrices, such that 'A = L * U'. With one\n",
+ " output argument Y, then the matrix returned by the LAPACK routines\n",
+ " is returned. If the input matrix is sparse then the matrix L is\n",
+ " embedded into U to give a return value similar to the full case.\n",
+ " For both full and sparse matrices, 'lu' loses the permutation\n",
+ " information.\n",
+ "\n",
+ " See also: luupdate, ilu, chol, hess, qr, qz, schur, svd.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help lu"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% time LU solution vs backslash\n",
+ "t_lu=zeros(100,1);\n",
+ "t_bs=zeros(100,1);\n",
+ "for N=1:100\n",
+ " A=rand(N,N);\n",
+ " y=rand(N,1);\n",
+ " [L,U,P]=lu(A);\n",
+ "\n",
+ " tic; d=inv(L)*y; x=inv(U)*d; t_lu(N)=toc;\n",
+ "\n",
+ " tic; x=inv(A)*y; t_bs(N)=toc;\n",
+ "end\n",
+ "plot([1:100],t_lu,[1:100],t_bs) \n",
+ "legend('LU decomp','Octave \\\\')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Spring-mass system\n",
+ "\n",
+ "4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "This matrix, K, is symmetric. \n",
+ "\n",
+ "`K(i,j)==K(j,i)`\n",
+ "\n",
+ "Now we can use,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "## 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": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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))"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/11_LU-decomposition/11_LU-decomposition.ipynb b/11_LU-decomposition/11_LU-decomposition.ipynb
new file mode 100644
index 0000000..726d345
--- /dev/null
+++ b/11_LU-decomposition/11_LU-decomposition.ipynb
@@ -0,0 +1,760 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# LU Decomposition\n",
+ "### efficient storage of matrices for solutions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Consider equation set:\n",
+ "\n",
+ "$y=Ax$\n",
+ "\n",
+ "Assume that we can perform Gauss elimination and achieve this formula:\n",
+ "\n",
+ "$Ux=d$ \n",
+ "\n",
+ "$U$: upper triangular matrix derived from Gauss elimination \n",
+ "\n",
+ "$d$: is the set of dependent variables after Gauss elimination. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Assume lower triangular matrix, $L$, exists \n",
+ "\n",
+ "- ones on the diagonal \n",
+ "- same dimensions as $U$ \n",
+ "- following is true:\n",
+ "\n",
+ "$L(Ux-d)=Ax-y=0$\n",
+ "\n",
+ "Now, $Ax=LUx$, so $A=LU$, and $y=Ld$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "$2x_{1}+x_{2}=1$\n",
+ "\n",
+ "$x_{1}+3x_{2}=1$\n",
+ "\n",
+ "\n",
+ "$\\left[ \\begin{array}{cc}\n",
+ "2 & 1 \\\\\n",
+ "1 & 3 \\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\end{array}\\right]=\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "1\\end{array}\\right]$\n",
+ "\n",
+ "f21=0.5\n",
+ "\n",
+ "A(2,1)=1-1 = 0 \n",
+ "\n",
+ "A(2,2)=3-0.5=2.5\n",
+ "\n",
+ "y(2)=1-0.5=0.5\n",
+ "\n",
+ "$L(Ux-d)=\n",
+ "\\left[ \\begin{array}{cc}\n",
+ "1 & 0 \\\\\n",
+ "0.5 & 1 \\end{array} \\right]\n",
+ "\\left(\\left[ \\begin{array}{cc}\n",
+ "2 & 1 \\\\\n",
+ "0 & 2.5 \\end{array} \\right]\n",
+ "\\left[\\begin{array}{c} \n",
+ "x_{1} \\\\ \n",
+ "x_{2} \\end{array}\\right]-\n",
+ "\\left[\\begin{array}{c} \n",
+ "1 \\\\\n",
+ "0.5\\end{array}\\right]\\right)=0$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "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": 13,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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": 14,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n",
+ "\n",
+ " -- Built-in Function: [L, U] = lu (A)\n",
+ " -- Built-in Function: [L, U, P] = lu (A)\n",
+ " -- Built-in Function: [L, U, P, Q] = lu (S)\n",
+ " -- Built-in Function: [L, U, P, Q, R] = lu (S)\n",
+ " -- Built-in Function: [...] = lu (S, THRES)\n",
+ " -- Built-in Function: Y = lu (...)\n",
+ " -- Built-in Function: [...] = lu (..., \"vector\")\n",
+ " Compute the LU decomposition of A.\n",
+ "\n",
+ " If A is full subroutines from LAPACK are used and if A is sparse\n",
+ " then UMFPACK is used.\n",
+ "\n",
+ " The result is returned in a permuted form, according to the\n",
+ " optional return value P. For example, given the matrix 'a = [1, 2;\n",
+ " 3, 4]',\n",
+ "\n",
+ " [l, u, p] = lu (A)\n",
+ "\n",
+ " returns\n",
+ "\n",
+ " l =\n",
+ "\n",
+ " 1.00000 0.00000\n",
+ " 0.33333 1.00000\n",
+ "\n",
+ " u =\n",
+ "\n",
+ " 3.00000 4.00000\n",
+ " 0.00000 0.66667\n",
+ "\n",
+ " p =\n",
+ "\n",
+ " 0 1\n",
+ " 1 0\n",
+ "\n",
+ " The matrix is not required to be square.\n",
+ "\n",
+ " When called with two or three output arguments and a spare input\n",
+ " matrix, 'lu' does not attempt to perform sparsity preserving column\n",
+ " permutations. Called with a fourth output argument, the sparsity\n",
+ " preserving column transformation Q is returned, such that 'P * A *\n",
+ " Q = L * U'.\n",
+ "\n",
+ " Called with a fifth output argument and a sparse input matrix, 'lu'\n",
+ " attempts to use a scaling factor R on the input matrix such that 'P\n",
+ " * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n",
+ " stable factorization.\n",
+ "\n",
+ " An additional input argument THRES, that defines the pivoting\n",
+ " threshold can be given. THRES can be a scalar, in which case it\n",
+ " defines the UMFPACK pivoting tolerance for both symmetric and\n",
+ " unsymmetric cases. If THRES is a 2-element vector, then the first\n",
+ " element defines the pivoting tolerance for the unsymmetric UMFPACK\n",
+ " pivoting strategy and the second for the symmetric strategy. By\n",
+ " default, the values defined by 'spparms' are used ([0.1, 0.001]).\n",
+ "\n",
+ " Given the string argument \"vector\", 'lu' returns the values of P\n",
+ " and Q as vector values, such that for full matrix, 'A (P,:) = L *\n",
+ " U', and 'R(P,:) * A (:, Q) = L * U'.\n",
+ "\n",
+ " With two output arguments, returns the permuted forms of the upper\n",
+ " and lower triangular matrices, such that 'A = L * U'. With one\n",
+ " output argument Y, then the matrix returned by the LAPACK routines\n",
+ " is returned. If the input matrix is sparse then the matrix L is\n",
+ " embedded into U to give a return value similar to the full case.\n",
+ " For both full and sparse matrices, 'lu' loses the permutation\n",
+ " information.\n",
+ "\n",
+ " See also: luupdate, ilu, chol, hess, qr, qz, schur, svd.\n",
+ "\n",
+ "Additional help for built-in functions and operators is\n",
+ "available in the online version of the manual. Use the command\n",
+ "'doc ' to search the manual index.\n",
+ "\n",
+ "Help and information about Octave is also available on the WWW\n",
+ "at http://www.octave.org and via the help@octave.org\n",
+ "mailing list.\n"
+ ]
+ }
+ ],
+ "source": [
+ "help lu"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% time LU solution vs backslash\n",
+ "t_lu=zeros(100,1);\n",
+ "t_bs=zeros(100,1);\n",
+ "for N=1:100\n",
+ " A=rand(N,N);\n",
+ " y=rand(N,1);\n",
+ " [L,U,P]=lu(A);\n",
+ "\n",
+ " tic; d=L\\y; x=U\\d; t_lu(N)=toc;\n",
+ "\n",
+ " tic; x=A\\y; t_bs(N)=toc;\n",
+ "end\n",
+ "plot([1:100],t_lu,[1:100],t_bs) \n",
+ "legend('LU decomp','Octave \\\\')\n",
+ "axis([0,100,0,0.001])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Spring-mass system\n",
+ "\n",
+ "4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
+ "\n",
+ "\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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "This matrix, K, is symmetric. \n",
+ "\n",
+ "`K(i,j)==K(j,i)`\n",
+ "\n",
+ "Now we can use,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "## 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": 12,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "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))"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Slideshow",
+ "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/11_LU-decomposition/LU_naive.m b/11_LU-decomposition/LU_naive.m
new file mode 100644
index 0000000..92efde6
--- /dev/null
+++ b/11_LU-decomposition/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/11_LU-decomposition/LU_pivot.m b/11_LU-decomposition/LU_pivot.m
new file mode 100644
index 0000000..37abb26
--- /dev/null
+++ b/11_LU-decomposition/LU_pivot.m
@@ -0,0 +1,36 @@
+function [L,U,P] = LU_pivot(A)
+% GaussPivot: Gauss elimination pivoting
+% x = GaussPivot(A,b): Gauss elimination with pivoting.
+% input:
+% A = coefficient matrix
+% b = 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;
+P=diag(ones(n,1));
+% forward elimination
+for k = 1:n-1
+ % partial pivoting
+ [big,i]=max(abs(U(k:n,k)));
+ ipr=i+k-1;
+ if ipr~=k
+ P([k,ipr],:)=P([ipr,k],:); % if the max is not the current index ipr, pivot count
+ %L([k,ipr],:)=L([ipr,k],:);
+ %U([k,ipr],:)=U([ipr,k],:);
+ end
+ 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/11_LU-decomposition/mass_springs.png b/11_LU-decomposition/mass_springs.png
new file mode 100644
index 0000000..259a333
Binary files /dev/null and b/11_LU-decomposition/mass_springs.png differ
diff --git a/11_LU-decomposition/mass_springs.svg b/11_LU-decomposition/mass_springs.svg
new file mode 100644
index 0000000..2bec8a4
--- /dev/null
+++ b/11_LU-decomposition/mass_springs.svg
@@ -0,0 +1,214 @@
+
+
+
+
diff --git a/11_LU-decomposition/octave-workspace b/11_LU-decomposition/octave-workspace
new file mode 100644
index 0000000..fdace26
Binary files /dev/null and b/11_LU-decomposition/octave-workspace differ