diff --git a/10_Gauss_elimination/4_dof.png b/10_Gauss_elimination/4_dof.png
new file mode 100644
index 0000000..caadeb5
Binary files /dev/null and b/10_Gauss_elimination/4_dof.png differ
diff --git a/10_Gauss_elimination/4_dof.svg b/10_Gauss_elimination/4_dof.svg
new file mode 100644
index 0000000..04bd409
--- /dev/null
+++ b/10_Gauss_elimination/4_dof.svg
@@ -0,0 +1,150 @@
+
+
+
+
diff --git a/10_Gauss_elimination/octave-workspace b/10_Gauss_elimination/octave-workspace
index 2a6ff70..8c437bb 100644
Binary files a/10_Gauss_elimination/octave-workspace and b/10_Gauss_elimination/octave-workspace differ
diff --git a/11_LU-decomposition/octave-workspace b/11_LU-decomposition/octave-workspace
index fdace26..8c437bb 100644
Binary files a/11_LU-decomposition/octave-workspace and b/11_LU-decomposition/octave-workspace differ
diff --git a/12_matrix_condition/.GS_rel.m.swp b/12_matrix_condition/.GS_rel.m.swp
new file mode 100644
index 0000000..58db6b3
Binary files /dev/null and b/12_matrix_condition/.GS_rel.m.swp differ
diff --git a/12_matrix_condition/.GaussSeidel.m.swp b/12_matrix_condition/.GaussSeidel.m.swp
new file mode 100644
index 0000000..43abe04
Binary files /dev/null and b/12_matrix_condition/.GaussSeidel.m.swp differ
diff --git a/12_matrix_condition/.Jacobi_rel.m.swp b/12_matrix_condition/.Jacobi_rel.m.swp
new file mode 100644
index 0000000..8f17e8f
Binary files /dev/null and b/12_matrix_condition/.Jacobi_rel.m.swp differ
diff --git a/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb b/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb
new file mode 100644
index 0000000..cc25231
--- /dev/null
+++ b/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb
@@ -0,0 +1,5261 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Condition of a matrix \n",
+ "### *just checked in to see what condition my condition was in*\n",
+ "### Matrix norms\n",
+ "\n",
+ "The Euclidean norm of a vector is measure of the magnitude (in 3D this would be: $|x|=\\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$) in general the equation is:\n",
+ "\n",
+ "$||x||_{e}=\\sqrt{\\sum_{i=1}^{n}x_{i}^{2}}$\n",
+ "\n",
+ "For a matrix, A, the same norm is called the Frobenius norm:\n",
+ "\n",
+ "$||A||_{f}=\\sqrt{\\sum_{i=1}^{n}\\sum_{j=1}^{m}A_{i,j}^{2}}$\n",
+ "\n",
+ "In general we can calculate any $p$-norm where\n",
+ "\n",
+ "$||A||_{p}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{p}}$\n",
+ "\n",
+ "so the p=1, 1-norm is \n",
+ "\n",
+ "$||A||_{1}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{1}}=\\sum_{i=1}^{n}\\sum_{i=1}^{m}|A_{i,j}|$\n",
+ "\n",
+ "$||A||_{\\infty}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{\\infty}}=\\max_{1\\le i \\le n}\\sum_{j=1}^{m}|A_{i,j}|$\n",
+ "\n",
+ "### Condition of Matrix\n",
+ "\n",
+ "The matrix condition is the product of \n",
+ "\n",
+ "$Cond(A) = ||A||\\cdot||A^{-1}||$ \n",
+ "\n",
+ "So each norm will have a different condition number, but the limit is $Cond(A)\\ge 1$\n",
+ "\n",
+ "An estimate of the rounding error is based on the condition of A:\n",
+ "\n",
+ "$\\frac{||\\Delta x||}{x} \\le Cond(A) \\frac{||\\Delta A||}{||A||}$\n",
+ "\n",
+ "So if the coefficients of A have accuracy to $10^{-t}\n",
+ "\n",
+ "and the condition of A, $Cond(A)=10^{c}$\n",
+ "\n",
+ "then the solution for x can have rounding errors up to $10^{c-t}$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 1.00000 0.50000 0.33333\n",
+ " 0.50000 0.33333 0.25000\n",
+ " 0.33333 0.25000 0.20000\n",
+ "\n",
+ "L =\n",
+ "\n",
+ " 1.00000 0.00000 0.00000\n",
+ " 0.50000 1.00000 0.00000\n",
+ " 0.33333 1.00000 1.00000\n",
+ "\n",
+ "U =\n",
+ "\n",
+ " 1.00000 0.50000 0.33333\n",
+ " 0.00000 0.08333 0.08333\n",
+ " 0.00000 -0.00000 0.00556\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5]\n",
+ "[L,U]=LU_naive(A)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "Then, $A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$\n",
+ "\n",
+ "$Ld_{1}=\\left[\\begin{array}{c}\n",
+ "1 \\\\\n",
+ "0 \\\\\n",
+ "0 \\end{array}\\right]$, $Ux_{1}=d_{1}$ ..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "invA =\n",
+ "\n",
+ " 9.0000 -36.0000 30.0000\n",
+ " -36.0000 192.0000 -180.0000\n",
+ " 30.0000 -180.0000 180.0000\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 1.0000e+00 3.5527e-15 2.9976e-15\n",
+ " -1.3249e-14 1.0000e+00 -9.1038e-15\n",
+ " 8.5117e-15 7.1054e-15 1.0000e+00\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "invA=zeros(3,3);\n",
+ "d1=L\\[1;0;0];\n",
+ "d2=L\\[0;1;0];\n",
+ "d3=L\\[0;0;1];\n",
+ "invA(:,1)=U\\d1; % shortcut invA(:,1)=A\\[1;0;0]\n",
+ "invA(:,2)=U\\d2;\n",
+ "invA(:,3)=U\\d3\n",
+ "invA*A"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Find the condition of A, $cond(A)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "normf_A = 1.4136\n",
+ "normf_invA = 372.21\n",
+ "cond_f_A = 526.16\n",
+ "ans = 1.4136\n",
+ "norm1_A = 1.8333\n",
+ "norm1_invA = 30.000\n",
+ "ans = 1.8333\n",
+ "cond_1_A = 55.000\n",
+ "norminf_A = 1.8333\n",
+ "norminf_invA = 30.000\n",
+ "ans = 1.8333\n",
+ "cond_inf_A = 55.000\n"
+ ]
+ }
+ ],
+ "source": [
+ "% Frobenius norm\n",
+ "normf_A = sqrt(sum(sum(A.^2)))\n",
+ "normf_invA = sqrt(sum(sum(invA.^2)))\n",
+ "\n",
+ "cond_f_A = normf_A*normf_invA\n",
+ "\n",
+ "norm(A,'fro')\n",
+ "\n",
+ "% p=1, column sum norm\n",
+ "norm1_A = max(sum(A,2))\n",
+ "norm1_invA = max(sum(invA,2))\n",
+ "norm(A,1)\n",
+ "\n",
+ "cond_1_A=norm1_A*norm1_invA\n",
+ "\n",
+ "% p=inf, row sum norm\n",
+ "norminf_A = max(sum(A,1))\n",
+ "norminf_invA = max(sum(invA,1))\n",
+ "norm(A,inf)\n",
+ "\n",
+ "cond_inf_A=norminf_A*norminf_invA\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with spring constants $K_{i}$. What does a high condition number mean for this problem? \n",
+ "\n",
+ "\n",
+ "\n",
+ "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n",
+ "\n",
+ "$m_{1}g+k_{2}(x_{2}-x_{1})-k_{1}x_{1}=0$\n",
+ "\n",
+ "$m_{2}g+k_{3}(x_{3}-x_{2})-k_{2}(x_{2}-x_{1})=0$\n",
+ "\n",
+ "$m_{3}g+k_{4}(x_{4}-x_{3})-k_{3}(x_{3}-x_{2})=0$\n",
+ "\n",
+ "$m_{4}g-k_{4}(x_{4}-x_{3})=0$\n",
+ "\n",
+ "in matrix form:\n",
+ "\n",
+ "$\\left[ \\begin{array}{cccc}\n",
+ "k_{1}+k_{2} & -k_{2} & 0 & 0 \\\\\n",
+ "-k_{2} & k_{2}+k_{3} & -k_{3} & 0 \\\\\n",
+ "0 & -k_{3} & k_{3}+k_{4} & -k_{4} \\\\\n",
+ "0 & 0 & -k_{4} & k_{4} \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\\\\n",
+ "x_{3} \\\\\n",
+ "x_{4} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "m_{1}g \\\\\n",
+ "m_{2}g \\\\\n",
+ "m_{3}g \\\\\n",
+ "m_{4}g \\end{array} \\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 100010 -100000 0 0\n",
+ " -100000 100010 -10 0\n",
+ " 0 -10 11 -1\n",
+ " 0 0 -1 1\n",
+ "\n",
+ "y =\n",
+ "\n",
+ " 9.8100\n",
+ " 19.6200\n",
+ " 29.4300\n",
+ " 39.2400\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "k1=10; % N/m\n",
+ "k2=100000;\n",
+ "k3=10;\n",
+ "k4=1;\n",
+ "m1=1; % kg\n",
+ "m2=2;\n",
+ "m3=3;\n",
+ "m4=4;\n",
+ "g=9.81; % m/s^2\n",
+ "K=[k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4]\n",
+ "y=[m1*g;m2*g;m3*g;m4*g]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 3.2004e+05\n",
+ "ans = 3.2004e+05\n",
+ "ans = 2.5925e+05\n",
+ "ans = 2.5293e+05\n"
+ ]
+ }
+ ],
+ "source": [
+ "cond(K,inf)\n",
+ "cond(K,1)\n",
+ "cond(K,'fro')\n",
+ "cond(K,2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "e =\n",
+ "\n",
+ " 7.9078e-01\n",
+ " 3.5881e+00\n",
+ " 1.7621e+01\n",
+ " 2.0001e+05\n",
+ "\n",
+ "ans = 2.5293e+05\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=eig(K)\n",
+ "max(e)/min(e)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## P=2 norm is ratio of biggest eigenvalue to smallest eigenvalue!\n",
+ "\n",
+ "no need to calculate the inv(K)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Iterative Methods\n",
+ "\n",
+ "## Gauss-Seidel method\n",
+ "\n",
+ "If we have an intial guess for each value of a vector $x$ that we are trying to solve, then it is easy enough to solve for one component given the others. \n",
+ "\n",
+ "Take a 3$\\times$3 matrix \n",
+ "\n",
+ "$Ax=b$\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc}\n",
+ "3 & -0.1 & -0.2 \\\\\n",
+ "0.1 & 7 & -0.3 \\\\\n",
+ "0.3 & -0.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",
+ "7.85 \\\\\n",
+ "-19.3 \\\\\n",
+ "71.4\\end{array} \\right]$\n",
+ "\n",
+ "$x_{1}=\\frac{7.85+0.1x_{2}+0.2x_{3}}{3}$\n",
+ "\n",
+ "$x_{2}=\\frac{-19.3-0.1x_{1}+0.3x_{3}}{7}$\n",
+ "\n",
+ "$x_{3}=\\frac{71.4+0.1x_{1}+0.2x_{2}}{10}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 3.00000 -0.10000 -0.20000\n",
+ " 0.10000 7.00000 -0.30000\n",
+ " 0.30000 -0.20000 10.00000\n",
+ "\n",
+ "b =\n",
+ "\n",
+ " 7.8500\n",
+ " -19.3000\n",
+ " 71.4000\n",
+ "\n",
+ "x =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[3 -0.1 -0.2;0.1 7 -0.3;0.3 -0.2 10]\n",
+ "b=[7.85;-19.3;71.4]\n",
+ "\n",
+ "x=A\\b"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "### Gauss-Seidel Iterative approach\n",
+ "\n",
+ "As a first guess, we can use $x_{1}=x_{2}=x_{3}=0$\n",
+ "\n",
+ "$x_{1}=\\frac{7.85+0.1(0)+0.3(0)}{3}=2.6167$\n",
+ "\n",
+ "$x_{2}=\\frac{-19.3-0.1(2.6167)+0.3(0)}{7}=-2.7945$\n",
+ "\n",
+ "$x_{3}=\\frac{71.4+0.1(2.6167)+0.2(-2.7945)}{10}=7.0056$\n",
+ "\n",
+ "Then, we update the guess:\n",
+ "\n",
+ "$x_{1}=\\frac{7.85+0.1(-2.7945)+0.3(7.0056)}{3}=2.9906$\n",
+ "\n",
+ "$x_{2}=\\frac{-19.3-0.1(2.9906)+0.3(7.0056)}{7}=-2.4996$\n",
+ "\n",
+ "$x_{3}=\\frac{71.4+0.1(2.9906)+0.2(-2.4966)}{10}=7.00029$\n",
+ "\n",
+ "The results are conveerging to the solution we found with `\\` of $x_{1}=3,~x_{2}=-2.5,~x_{3}=7$\n",
+ "\n",
+ "We could also use an iterative method that solves for all of the x-components in one step:\n",
+ "\n",
+ "### Jacobi method\n",
+ "\n",
+ "$x_{1}^{i}=\\frac{7.85+0.1x_{2}^{i-1}+0.3x_{3}^{i-1}}{3}$\n",
+ "\n",
+ "$x_{2}^{i}=\\frac{-19.3-0.1x_{1}^{i-1}+0.3x_{3}^{i-1}}{7}$\n",
+ "\n",
+ "$x_{3}^{i}=\\frac{71.4+0.1x_{1}^{i-1}+0.2x_{2}^{i-1}}{10}$\n",
+ "\n",
+ "Here the solution is a matrix multiplication and vector addition\n",
+ "\n",
+ "$\\left[ \\begin{array}{c}\n",
+ "x_{1}^{i} \\\\\n",
+ "x_{2}^{i} \\\\\n",
+ "x_{3}^{i} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "7.85/3 \\\\\n",
+ "-19.3/7 \\\\\n",
+ "71.4/10\\end{array} \\right]-\n",
+ "\\left[ \\begin{array}{ccc}\n",
+ "0 & 0.1/3 & 0.2/3 \\\\\n",
+ "0.1/7 & 0 & -0.3/7 \\\\\n",
+ "0.3/10 & -0.2/10 & 0 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1}^{i-1} \\\\\n",
+ "x_{2}^{i-1} \\\\\n",
+ "x_{3}^{i-1} \\end{array} \\right]$\n",
+ "\n",
+ "|x_{j}|Jacobi method |vs| Gauss-Seidel |\n",
+ "|--------|------------------------------|---|-------------------------------|\n",
+ "|$x_{1}^{i}=$ | $\\frac{7.85+0.1x_{2}^{i-1}+0.3x_{3}^{i-1}}{3}$ | | $\\frac{7.85+0.1x_{2}^{i-1}+0.3x_{3}^{i-1}}{3}$|\n",
+ "|$x_{2}^{i}=$ | $\\frac{-19.3-0.1x_{1}^{i-1}+0.3x_{3}^{i-1}}{7}$ | | $\\frac{-19.3-0.1x_{1}^{i}+0.3x_{3}^{i-1}}{7}$ |\n",
+ "|$x_{3}^{i}=$ | $\\frac{71.4+0.1x_{1}^{i-1}+0.2x_{2}^{i-1}}{10}$ | | $\\frac{71.4+0.1x_{1}^{i}+0.2x_{2}^{i}}{10}$|"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ba =\n",
+ "\n",
+ " 2.6167\n",
+ " -2.7571\n",
+ " 7.1400\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0.00000 -0.10000 -0.20000\n",
+ " 0.10000 0.00000 -0.30000\n",
+ " 0.30000 -0.20000 0.00000\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0.000000 -0.033333 -0.066667\n",
+ " 0.014286 0.000000 -0.042857\n",
+ " 0.030000 -0.020000 0.000000\n",
+ "\n",
+ "x1 =\n",
+ "\n",
+ " 2.6167\n",
+ " -2.7571\n",
+ " 7.1400\n",
+ "\n",
+ "x2 =\n",
+ "\n",
+ " 3.0008\n",
+ " -2.4885\n",
+ " 7.0064\n",
+ "\n",
+ "x3 =\n",
+ "\n",
+ " 3.0008\n",
+ " -2.4997\n",
+ " 7.0002\n",
+ "\n",
+ "solution is converging to [3,-2.5,7]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "ba=b./diag(A) % or ba=b./[A(1,1);A(2,2);A(3,3)]\n",
+ "sA=A-diag(diag(A)) % A with zeros on diagonal\n",
+ "sA(1,:)=sA(1,:)/A(1,1);\n",
+ "sA(2,:)=sA(2,:)/A(2,2);\n",
+ "sA(3,:)=sA(3,:)/A(3,3)\n",
+ "x0=[0;0;0];\n",
+ "x1=ba-sA*x0\n",
+ "x2=ba-sA*x1\n",
+ "x3=ba-sA*x2\n",
+ "fprintf('solution is converging to [3,-2.5,7]]\\n')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 3\n",
+ " 7\n",
+ " 10\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ "Diagonal Matrix\n",
+ "\n",
+ " 3 0 0\n",
+ " 0 7 0\n",
+ " 0 0 10\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "diag(A)\n",
+ "diag(diag(A))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This method works if problem is diagonally dominant, \n",
+ "\n",
+ "$|a_{ii}|>\\sum_{j=1,j\\ne i}^{n}|a_{ij}|$\n",
+ "\n",
+ "If this condition is true, then Jacobi or Gauss-Seidel should converge\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 0.10000 1.00000 3.00000\n",
+ " 1.00000 0.20000 3.00000\n",
+ " 5.00000 2.00000 0.30000\n",
+ "\n",
+ "b =\n",
+ "\n",
+ " 12\n",
+ " 2\n",
+ " 4\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -2.9393\n",
+ " 9.1933\n",
+ " 1.0336\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[0.1,1,3;1,0.2,3;5,2,0.3]\n",
+ "b=[12;2;4]\n",
+ "A\\b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ba =\n",
+ "\n",
+ " 120.000\n",
+ " 10.000\n",
+ " 13.333\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0 1 3\n",
+ " 1 0 3\n",
+ " 5 2 0\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0.00000 10.00000 30.00000\n",
+ " 5.00000 0.00000 15.00000\n",
+ " 16.66667 6.66667 0.00000\n",
+ "\n",
+ "x1 =\n",
+ "\n",
+ " 120.000\n",
+ " 10.000\n",
+ " 13.333\n",
+ "\n",
+ "x2 =\n",
+ "\n",
+ " -380.00\n",
+ " -790.00\n",
+ " -2053.33\n",
+ "\n",
+ "x3 =\n",
+ "\n",
+ " 6.9620e+04\n",
+ " 3.2710e+04\n",
+ " 1.1613e+04\n",
+ "\n",
+ "solution is not converging to [-2.93,9.19,1.03]\n"
+ ]
+ }
+ ],
+ "source": [
+ "ba=b./diag(A) % or ba=b./[A(1,1);A(2,2);A(3,3)]\n",
+ "sA=A-diag(diag(A)) % A with zeros on diagonal\n",
+ "sA(1,:)=sA(1,:)/A(1,1);\n",
+ "sA(2,:)=sA(2,:)/A(2,2);\n",
+ "sA(3,:)=sA(3,:)/A(3,3)\n",
+ "x0=[0;0;0];\n",
+ "x1=ba-sA*x0\n",
+ "x2=ba-sA*x1\n",
+ "x3=ba-sA*x2\n",
+ "fprintf('solution is not converging to [-2.93,9.19,1.03]\\n')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Gauss-Seidel with Relaxation\n",
+ "\n",
+ "In order to force the solution to converge faster, we can introduce a relaxation term $\\lambda$. \n",
+ "\n",
+ "where the new x values are weighted between the old and new:\n",
+ "\n",
+ "$x^{i}=\\lambda x^{i}+(1-\\lambda)x^{i-1}$\n",
+ "\n",
+ "after solving for x, lambda weights the current approximation with the previous approximation for the updated x\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 3.00000 -0.10000 -0.20000\n",
+ " 0.10000 7.00000 -0.30000\n",
+ " 0.30000 -0.20000 10.00000\n",
+ "\n",
+ "b =\n",
+ "\n",
+ " 7.8500\n",
+ " -19.3000\n",
+ " 71.4000\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% rearrange A and b\n",
+ "A=[3 -0.1 -0.2;0.1 7 -0.3;0.3 -0.2 10]\n",
+ "b=[7.85;-19.3;71.4]\n",
+ "\n",
+ "iters=zeros(100,1);\n",
+ "for i=1:100\n",
+ " lambda=2/100*i;\n",
+ " [x,ea,iters(i)]=Jacobi_rel(A,b,lambda);\n",
+ "end\n",
+ "plot([1:100]*2/100,iters) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 107,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "l = 0.99158\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "l=fminbnd(@(l) lambda_fcn(A,b,l),0.5,1.5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 108,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A\\b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 109,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 1.8289e-07\n",
+ " 2.1984e-08\n",
+ " 2.3864e-08\n",
+ "\n",
+ "iter = 8\n",
+ "x =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 1.9130e-08\n",
+ " 7.6449e-08\n",
+ " 3.3378e-08\n",
+ "\n",
+ "iter = 8\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,ea,iter]=Jacobi_rel(A,b,l,0.000001)\n",
+ "[x,ea,iter]=Jacobi_rel(A,b,1,0.000001)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Nonlinear Systems\n",
+ "\n",
+ "Consider two simultaneous nonlinear equations with two unknowns:\n",
+ "\n",
+ "$x_{1}^{2}+x_{1}x_{2}=10$\n",
+ "\n",
+ "$x_{2}+3x_{1}x_{2}^{2}=57$\n",
+ "\n",
+ "Graphically, we are looking for the solution:\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x11=linspace(0.5,3);\n",
+ "x12=(10-x11.^2)./x11;\n",
+ "\n",
+ "x22=linspace(2,8);\n",
+ "x21=(57-x22).*x22.^-2/3;\n",
+ "\n",
+ "plot(x11,x12,x21,x22)\n",
+ "% Solution at x_1=2, x_2=3\n",
+ "hold on;\n",
+ "plot(2,3,'o')\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Newton-Raphson part II\n",
+ "\n",
+ "Remember the first order approximation for the next point in a function is:\n",
+ "\n",
+ "$f(x_{i+1})=f(x_{i})+(x_{i+1}-x_{i})f'(x_{i})$\n",
+ "\n",
+ "then, $f(x_{i+1})=0$ so we are left with:\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n",
+ "\n",
+ "We can use the same formula, but now we have multiple dimensions so we need to determine the Jacobian\n",
+ "\n",
+ "$[J]=\\left[ \\begin{array}{cccc}\n",
+ "\\frac{\\partial f_{1,i}}{\\partial x_{1}} & \\frac{\\partial f_{1,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{1,i}}{\\partial x_{n}} \\\\\n",
+ "\\frac{\\partial f_{2,i}}{\\partial x_{1}} & \\frac{\\partial f_{2,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{2,i}}{\\partial x_{n}} \\\\\n",
+ "\\vdots & \\vdots & & \\vdots \\\\\n",
+ "\\frac{\\partial f_{n,i}}{\\partial x_{1}} & \\frac{\\partial f_{n,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{n,i}}{\\partial x_{n}} \\\\\n",
+ "\\end{array} \\right]$\n",
+ "\n",
+ "$\\left[ \\begin{array}{c}\n",
+ "f_{1,i+1} \\\\\n",
+ "f_{2,i+1} \\\\\n",
+ "\\vdots \\\\\n",
+ "f_{n,i+1}\\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "f_{1,i} \\\\\n",
+ "f_{2,i} \\\\\n",
+ "\\vdots \\\\\n",
+ "f_{n,i}\\end{array} \\right]+\n",
+ "\\left[ \\begin{array}{cccc}\n",
+ "\\frac{\\partial f_{1,i}}{\\partial x_{1}} & \\frac{\\partial f_{1,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{1,i}}{\\partial x_{n}} \\\\\n",
+ "\\frac{\\partial f_{2,i}}{\\partial x_{1}} & \\frac{\\partial f_{2,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{2,i}}{\\partial x_{n}} \\\\\n",
+ "\\vdots & \\vdots & & \\vdots \\\\\n",
+ "\\frac{\\partial f_{n,i}}{\\partial x_{1}} & \\frac{\\partial f_{n,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{n,i}}{\\partial x_{n}} \\\\\n",
+ "\\end{array} \\right]\n",
+ "\\left( \\left[ \\begin{array}{c}\n",
+ "x_{i+1} \\\\\n",
+ "x_{i+1} \\\\\n",
+ "\\vdots \\\\\n",
+ "x_{i+1}\\end{array} \\right]-\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1,i} \\\\\n",
+ "x_{2,i} \\\\\n",
+ "\\vdots \\\\\n",
+ "x_{n,i}\\end{array} \\right]\\right)$\n",
+ "\n",
+ "### Solution is again in the form Ax=b\n",
+ "\n",
+ "$[J]([x_{i+1}]-[x_{i}])=-[f]$\n",
+ "\n",
+ "so\n",
+ "\n",
+ "$[x_{i+1}]= [x_{i}]-[J]^{-1}[f]$\n",
+ "\n",
+ "## Example of Jacobian calculation\n",
+ "\n",
+ "### Nonlinear springs supporting two masses in series\n",
+ "\n",
+ "Two springs are connected to two masses, with $m_1$=1 kg and $m_{2}$=2 kg. The springs are identical, but they have nonlinear spring constants, of $k_1$=100 N/m and $k_2$=-10 N/m\n",
+ "\n",
+ "We want to solve for the final position of the masses ($x_1$ and $x_2$)\n",
+ "\n",
+ "$m_{1}g+k_{1}(x_{2}-x_{1})+k_{2}(x_{2}-x_{1})^{2}+k_{1}x_{1}+k_{2}x_{1}^{2}=0$\n",
+ "\n",
+ "$m_{2}g-k_{1}(x_{2}-x_{1})-k_{2}(x_2-x_1)^{2}=0$\n",
+ "\n",
+ "$J(1,1)=\\frac{\\partial f_{1}}{\\partial x_{1}}=-k_{1}-2k_{2}(x_{2}-x_{1})+k_{1}+2k_{2}x_{1}$\n",
+ "\n",
+ "$J(1,2)=\\frac{\\partial f_1}{\\partial x_{2}}=k_{1}+2k_{2}(x_{2}-x_{1})$\n",
+ "\n",
+ "$J(2,1)=\\frac{\\partial f_2}{\\partial x_{1}}=k_{1}+2k_{2}(x_{2}-x_{1})$\n",
+ "\n",
+ "$J(2,2)=\\frac{\\partial f_2}{\\partial x_{2}}=-k_{1}-2k_{2}(x_{2}-x_{1})$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "m1=1; % kg \n",
+ "m2=2; % kg\n",
+ "k1=100; % N/m\n",
+ "k2=-10; % N/m^2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "function [f,J]=mass_spring(x)\n",
+ " % Function to calculate function values f1 and f2 as well as Jacobian \n",
+ " % for 2 masses and 2 identical nonlinear springs\n",
+ " m1=1; % kg \n",
+ " m2=2; % kg\n",
+ " k1=100; % N/m\n",
+ " k2=-10; % N/m^2\n",
+ " g=9.81; % m/s^2\n",
+ " x1=x(1);\n",
+ " x2=x(2);\n",
+ " J=[-k1-2*k2*(x2-x1)-k1-2*k2*x1,k1+2*k2*(x2-x1);\n",
+ " k1+2*k2*(x2-x1),-k1-2*k2*(x2-x1)];\n",
+ " f=[m1*g+k1*(x2-x1)+k2*(x2-x1).^2-k1*x1-k2*x1^2;\n",
+ " m2*g-k1*(x2-x1)-k2*(x2-x1).^2];\n",
+ "end\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "f =\n",
+ "\n",
+ " -190.19\n",
+ " 129.62\n",
+ "\n",
+ "J =\n",
+ "\n",
+ " -200 120\n",
+ " 120 -120\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "[f,J]=mass_spring([1,0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x1 =\n",
+ "\n",
+ " -1.5142\n",
+ " -1.4341\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 2.9812\n",
+ " 2.3946\n",
+ "\n",
+ "x2 =\n",
+ "\n",
+ " 0.049894\n",
+ " 0.248638\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 31.3492\n",
+ " 6.7678\n",
+ "\n",
+ "x3 =\n",
+ "\n",
+ " 0.29701\n",
+ " 0.49722\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 0.83201\n",
+ " 0.49995\n",
+ "\n",
+ "x =\n",
+ "\n",
+ " 0.29701\n",
+ " 0.49722\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 0.021392\n",
+ " 0.012890\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 1.4786e-05\n",
+ " 8.9091e-06\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 7.0642e-12\n",
+ " 4.2565e-12\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x0=[3;2];\n",
+ "[f0,J0]=mass_spring(x0);\n",
+ "x1=x0-J0\\f0\n",
+ "ea=(x1-x0)./x1\n",
+ "[f1,J1]=mass_spring(x1);\n",
+ "x2=x1-J1\\f1\n",
+ "ea=(x2-x1)./x2\n",
+ "[f2,J2]=mass_spring(x2);\n",
+ "x3=x2-J2\\f2\n",
+ "ea=(x3-x2)./x3\n",
+ "x=x3\n",
+ "for i=1:3\n",
+ " xold=x;\n",
+ " [f,J]=mass_spring(x);\n",
+ " x=x-J\\f;\n",
+ " ea=(x-xold)./x\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.30351\n",
+ " 0.50372\n",
+ "\n",
+ "X0 =\n",
+ "\n",
+ " 0.30351\n",
+ " 0.50372\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x\n",
+ "X0=fsolve(@(x) mass_spring(x),[3;5])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[X,Y]=meshgrid(linspace(0,10,20),linspace(0,10,20));\n",
+ "[N,M]=size(X);\n",
+ "F=zeros(size(X));\n",
+ "for i=1:N\n",
+ " for j=1:M\n",
+ " [f,~]=mass_spring([X(i,j),Y(i,j)]);\n",
+ " F1(i,j)=f(1);\n",
+ " F2(i,j)=f(2);\n",
+ " end\n",
+ "end\n",
+ "mesh(X,Y,F1)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')\n",
+ "colorbar()\n",
+ "figure()\n",
+ "mesh(X,Y,F2)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')\n",
+ "colorbar()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/12_matrix_condition/.lambda_fcn.m.swp b/12_matrix_condition/.lambda_fcn.m.swp
new file mode 100644
index 0000000..5031365
Binary files /dev/null and b/12_matrix_condition/.lambda_fcn.m.swp differ
diff --git a/12_matrix_condition/12_matrix_condition.ipynb b/12_matrix_condition/12_matrix_condition.ipynb
new file mode 100644
index 0000000..cc25231
--- /dev/null
+++ b/12_matrix_condition/12_matrix_condition.ipynb
@@ -0,0 +1,5261 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Condition of a matrix \n",
+ "### *just checked in to see what condition my condition was in*\n",
+ "### Matrix norms\n",
+ "\n",
+ "The Euclidean norm of a vector is measure of the magnitude (in 3D this would be: $|x|=\\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$) in general the equation is:\n",
+ "\n",
+ "$||x||_{e}=\\sqrt{\\sum_{i=1}^{n}x_{i}^{2}}$\n",
+ "\n",
+ "For a matrix, A, the same norm is called the Frobenius norm:\n",
+ "\n",
+ "$||A||_{f}=\\sqrt{\\sum_{i=1}^{n}\\sum_{j=1}^{m}A_{i,j}^{2}}$\n",
+ "\n",
+ "In general we can calculate any $p$-norm where\n",
+ "\n",
+ "$||A||_{p}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{p}}$\n",
+ "\n",
+ "so the p=1, 1-norm is \n",
+ "\n",
+ "$||A||_{1}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{1}}=\\sum_{i=1}^{n}\\sum_{i=1}^{m}|A_{i,j}|$\n",
+ "\n",
+ "$||A||_{\\infty}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{\\infty}}=\\max_{1\\le i \\le n}\\sum_{j=1}^{m}|A_{i,j}|$\n",
+ "\n",
+ "### Condition of Matrix\n",
+ "\n",
+ "The matrix condition is the product of \n",
+ "\n",
+ "$Cond(A) = ||A||\\cdot||A^{-1}||$ \n",
+ "\n",
+ "So each norm will have a different condition number, but the limit is $Cond(A)\\ge 1$\n",
+ "\n",
+ "An estimate of the rounding error is based on the condition of A:\n",
+ "\n",
+ "$\\frac{||\\Delta x||}{x} \\le Cond(A) \\frac{||\\Delta A||}{||A||}$\n",
+ "\n",
+ "So if the coefficients of A have accuracy to $10^{-t}\n",
+ "\n",
+ "and the condition of A, $Cond(A)=10^{c}$\n",
+ "\n",
+ "then the solution for x can have rounding errors up to $10^{c-t}$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 1.00000 0.50000 0.33333\n",
+ " 0.50000 0.33333 0.25000\n",
+ " 0.33333 0.25000 0.20000\n",
+ "\n",
+ "L =\n",
+ "\n",
+ " 1.00000 0.00000 0.00000\n",
+ " 0.50000 1.00000 0.00000\n",
+ " 0.33333 1.00000 1.00000\n",
+ "\n",
+ "U =\n",
+ "\n",
+ " 1.00000 0.50000 0.33333\n",
+ " 0.00000 0.08333 0.08333\n",
+ " 0.00000 -0.00000 0.00556\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5]\n",
+ "[L,U]=LU_naive(A)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "Then, $A^{-1}=(LU)^{-1}=U^{-1}L^{-1}$\n",
+ "\n",
+ "$Ld_{1}=\\left[\\begin{array}{c}\n",
+ "1 \\\\\n",
+ "0 \\\\\n",
+ "0 \\end{array}\\right]$, $Ux_{1}=d_{1}$ ..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "invA =\n",
+ "\n",
+ " 9.0000 -36.0000 30.0000\n",
+ " -36.0000 192.0000 -180.0000\n",
+ " 30.0000 -180.0000 180.0000\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 1.0000e+00 3.5527e-15 2.9976e-15\n",
+ " -1.3249e-14 1.0000e+00 -9.1038e-15\n",
+ " 8.5117e-15 7.1054e-15 1.0000e+00\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "invA=zeros(3,3);\n",
+ "d1=L\\[1;0;0];\n",
+ "d2=L\\[0;1;0];\n",
+ "d3=L\\[0;0;1];\n",
+ "invA(:,1)=U\\d1; % shortcut invA(:,1)=A\\[1;0;0]\n",
+ "invA(:,2)=U\\d2;\n",
+ "invA(:,3)=U\\d3\n",
+ "invA*A"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Find the condition of A, $cond(A)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "normf_A = 1.4136\n",
+ "normf_invA = 372.21\n",
+ "cond_f_A = 526.16\n",
+ "ans = 1.4136\n",
+ "norm1_A = 1.8333\n",
+ "norm1_invA = 30.000\n",
+ "ans = 1.8333\n",
+ "cond_1_A = 55.000\n",
+ "norminf_A = 1.8333\n",
+ "norminf_invA = 30.000\n",
+ "ans = 1.8333\n",
+ "cond_inf_A = 55.000\n"
+ ]
+ }
+ ],
+ "source": [
+ "% Frobenius norm\n",
+ "normf_A = sqrt(sum(sum(A.^2)))\n",
+ "normf_invA = sqrt(sum(sum(invA.^2)))\n",
+ "\n",
+ "cond_f_A = normf_A*normf_invA\n",
+ "\n",
+ "norm(A,'fro')\n",
+ "\n",
+ "% p=1, column sum norm\n",
+ "norm1_A = max(sum(A,2))\n",
+ "norm1_invA = max(sum(invA,2))\n",
+ "norm(A,1)\n",
+ "\n",
+ "cond_1_A=norm1_A*norm1_invA\n",
+ "\n",
+ "% p=inf, row sum norm\n",
+ "norminf_A = max(sum(A,1))\n",
+ "norminf_invA = max(sum(invA,1))\n",
+ "norm(A,inf)\n",
+ "\n",
+ "cond_inf_A=norminf_A*norminf_invA\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with spring constants $K_{i}$. What does a high condition number mean for this problem? \n",
+ "\n",
+ "\n",
+ "\n",
+ "The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n",
+ "\n",
+ "$m_{1}g+k_{2}(x_{2}-x_{1})-k_{1}x_{1}=0$\n",
+ "\n",
+ "$m_{2}g+k_{3}(x_{3}-x_{2})-k_{2}(x_{2}-x_{1})=0$\n",
+ "\n",
+ "$m_{3}g+k_{4}(x_{4}-x_{3})-k_{3}(x_{3}-x_{2})=0$\n",
+ "\n",
+ "$m_{4}g-k_{4}(x_{4}-x_{3})=0$\n",
+ "\n",
+ "in matrix form:\n",
+ "\n",
+ "$\\left[ \\begin{array}{cccc}\n",
+ "k_{1}+k_{2} & -k_{2} & 0 & 0 \\\\\n",
+ "-k_{2} & k_{2}+k_{3} & -k_{3} & 0 \\\\\n",
+ "0 & -k_{3} & k_{3}+k_{4} & -k_{4} \\\\\n",
+ "0 & 0 & -k_{4} & k_{4} \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1} \\\\\n",
+ "x_{2} \\\\\n",
+ "x_{3} \\\\\n",
+ "x_{4} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "m_{1}g \\\\\n",
+ "m_{2}g \\\\\n",
+ "m_{3}g \\\\\n",
+ "m_{4}g \\end{array} \\right]$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 100010 -100000 0 0\n",
+ " -100000 100010 -10 0\n",
+ " 0 -10 11 -1\n",
+ " 0 0 -1 1\n",
+ "\n",
+ "y =\n",
+ "\n",
+ " 9.8100\n",
+ " 19.6200\n",
+ " 29.4300\n",
+ " 39.2400\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "k1=10; % N/m\n",
+ "k2=100000;\n",
+ "k3=10;\n",
+ "k4=1;\n",
+ "m1=1; % kg\n",
+ "m2=2;\n",
+ "m3=3;\n",
+ "m4=4;\n",
+ "g=9.81; % m/s^2\n",
+ "K=[k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4]\n",
+ "y=[m1*g;m2*g;m3*g;m4*g]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 3.2004e+05\n",
+ "ans = 3.2004e+05\n",
+ "ans = 2.5925e+05\n",
+ "ans = 2.5293e+05\n"
+ ]
+ }
+ ],
+ "source": [
+ "cond(K,inf)\n",
+ "cond(K,1)\n",
+ "cond(K,'fro')\n",
+ "cond(K,2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "e =\n",
+ "\n",
+ " 7.9078e-01\n",
+ " 3.5881e+00\n",
+ " 1.7621e+01\n",
+ " 2.0001e+05\n",
+ "\n",
+ "ans = 2.5293e+05\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=eig(K)\n",
+ "max(e)/min(e)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## P=2 norm is ratio of biggest eigenvalue to smallest eigenvalue!\n",
+ "\n",
+ "no need to calculate the inv(K)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Iterative Methods\n",
+ "\n",
+ "## Gauss-Seidel method\n",
+ "\n",
+ "If we have an intial guess for each value of a vector $x$ that we are trying to solve, then it is easy enough to solve for one component given the others. \n",
+ "\n",
+ "Take a 3$\\times$3 matrix \n",
+ "\n",
+ "$Ax=b$\n",
+ "\n",
+ "$\\left[ \\begin{array}{ccc}\n",
+ "3 & -0.1 & -0.2 \\\\\n",
+ "0.1 & 7 & -0.3 \\\\\n",
+ "0.3 & -0.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",
+ "7.85 \\\\\n",
+ "-19.3 \\\\\n",
+ "71.4\\end{array} \\right]$\n",
+ "\n",
+ "$x_{1}=\\frac{7.85+0.1x_{2}+0.2x_{3}}{3}$\n",
+ "\n",
+ "$x_{2}=\\frac{-19.3-0.1x_{1}+0.3x_{3}}{7}$\n",
+ "\n",
+ "$x_{3}=\\frac{71.4+0.1x_{1}+0.2x_{2}}{10}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 3.00000 -0.10000 -0.20000\n",
+ " 0.10000 7.00000 -0.30000\n",
+ " 0.30000 -0.20000 10.00000\n",
+ "\n",
+ "b =\n",
+ "\n",
+ " 7.8500\n",
+ " -19.3000\n",
+ " 71.4000\n",
+ "\n",
+ "x =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[3 -0.1 -0.2;0.1 7 -0.3;0.3 -0.2 10]\n",
+ "b=[7.85;-19.3;71.4]\n",
+ "\n",
+ "x=A\\b"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "### Gauss-Seidel Iterative approach\n",
+ "\n",
+ "As a first guess, we can use $x_{1}=x_{2}=x_{3}=0$\n",
+ "\n",
+ "$x_{1}=\\frac{7.85+0.1(0)+0.3(0)}{3}=2.6167$\n",
+ "\n",
+ "$x_{2}=\\frac{-19.3-0.1(2.6167)+0.3(0)}{7}=-2.7945$\n",
+ "\n",
+ "$x_{3}=\\frac{71.4+0.1(2.6167)+0.2(-2.7945)}{10}=7.0056$\n",
+ "\n",
+ "Then, we update the guess:\n",
+ "\n",
+ "$x_{1}=\\frac{7.85+0.1(-2.7945)+0.3(7.0056)}{3}=2.9906$\n",
+ "\n",
+ "$x_{2}=\\frac{-19.3-0.1(2.9906)+0.3(7.0056)}{7}=-2.4996$\n",
+ "\n",
+ "$x_{3}=\\frac{71.4+0.1(2.9906)+0.2(-2.4966)}{10}=7.00029$\n",
+ "\n",
+ "The results are conveerging to the solution we found with `\\` of $x_{1}=3,~x_{2}=-2.5,~x_{3}=7$\n",
+ "\n",
+ "We could also use an iterative method that solves for all of the x-components in one step:\n",
+ "\n",
+ "### Jacobi method\n",
+ "\n",
+ "$x_{1}^{i}=\\frac{7.85+0.1x_{2}^{i-1}+0.3x_{3}^{i-1}}{3}$\n",
+ "\n",
+ "$x_{2}^{i}=\\frac{-19.3-0.1x_{1}^{i-1}+0.3x_{3}^{i-1}}{7}$\n",
+ "\n",
+ "$x_{3}^{i}=\\frac{71.4+0.1x_{1}^{i-1}+0.2x_{2}^{i-1}}{10}$\n",
+ "\n",
+ "Here the solution is a matrix multiplication and vector addition\n",
+ "\n",
+ "$\\left[ \\begin{array}{c}\n",
+ "x_{1}^{i} \\\\\n",
+ "x_{2}^{i} \\\\\n",
+ "x_{3}^{i} \\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "7.85/3 \\\\\n",
+ "-19.3/7 \\\\\n",
+ "71.4/10\\end{array} \\right]-\n",
+ "\\left[ \\begin{array}{ccc}\n",
+ "0 & 0.1/3 & 0.2/3 \\\\\n",
+ "0.1/7 & 0 & -0.3/7 \\\\\n",
+ "0.3/10 & -0.2/10 & 0 \\end{array} \\right]\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1}^{i-1} \\\\\n",
+ "x_{2}^{i-1} \\\\\n",
+ "x_{3}^{i-1} \\end{array} \\right]$\n",
+ "\n",
+ "|x_{j}|Jacobi method |vs| Gauss-Seidel |\n",
+ "|--------|------------------------------|---|-------------------------------|\n",
+ "|$x_{1}^{i}=$ | $\\frac{7.85+0.1x_{2}^{i-1}+0.3x_{3}^{i-1}}{3}$ | | $\\frac{7.85+0.1x_{2}^{i-1}+0.3x_{3}^{i-1}}{3}$|\n",
+ "|$x_{2}^{i}=$ | $\\frac{-19.3-0.1x_{1}^{i-1}+0.3x_{3}^{i-1}}{7}$ | | $\\frac{-19.3-0.1x_{1}^{i}+0.3x_{3}^{i-1}}{7}$ |\n",
+ "|$x_{3}^{i}=$ | $\\frac{71.4+0.1x_{1}^{i-1}+0.2x_{2}^{i-1}}{10}$ | | $\\frac{71.4+0.1x_{1}^{i}+0.2x_{2}^{i}}{10}$|"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ba =\n",
+ "\n",
+ " 2.6167\n",
+ " -2.7571\n",
+ " 7.1400\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0.00000 -0.10000 -0.20000\n",
+ " 0.10000 0.00000 -0.30000\n",
+ " 0.30000 -0.20000 0.00000\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0.000000 -0.033333 -0.066667\n",
+ " 0.014286 0.000000 -0.042857\n",
+ " 0.030000 -0.020000 0.000000\n",
+ "\n",
+ "x1 =\n",
+ "\n",
+ " 2.6167\n",
+ " -2.7571\n",
+ " 7.1400\n",
+ "\n",
+ "x2 =\n",
+ "\n",
+ " 3.0008\n",
+ " -2.4885\n",
+ " 7.0064\n",
+ "\n",
+ "x3 =\n",
+ "\n",
+ " 3.0008\n",
+ " -2.4997\n",
+ " 7.0002\n",
+ "\n",
+ "solution is converging to [3,-2.5,7]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "ba=b./diag(A) % or ba=b./[A(1,1);A(2,2);A(3,3)]\n",
+ "sA=A-diag(diag(A)) % A with zeros on diagonal\n",
+ "sA(1,:)=sA(1,:)/A(1,1);\n",
+ "sA(2,:)=sA(2,:)/A(2,2);\n",
+ "sA(3,:)=sA(3,:)/A(3,3)\n",
+ "x0=[0;0;0];\n",
+ "x1=ba-sA*x0\n",
+ "x2=ba-sA*x1\n",
+ "x3=ba-sA*x2\n",
+ "fprintf('solution is converging to [3,-2.5,7]]\\n')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 3\n",
+ " 7\n",
+ " 10\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ "Diagonal Matrix\n",
+ "\n",
+ " 3 0 0\n",
+ " 0 7 0\n",
+ " 0 0 10\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "diag(A)\n",
+ "diag(diag(A))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This method works if problem is diagonally dominant, \n",
+ "\n",
+ "$|a_{ii}|>\\sum_{j=1,j\\ne i}^{n}|a_{ij}|$\n",
+ "\n",
+ "If this condition is true, then Jacobi or Gauss-Seidel should converge\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 0.10000 1.00000 3.00000\n",
+ " 1.00000 0.20000 3.00000\n",
+ " 5.00000 2.00000 0.30000\n",
+ "\n",
+ "b =\n",
+ "\n",
+ " 12\n",
+ " 2\n",
+ " 4\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -2.9393\n",
+ " 9.1933\n",
+ " 1.0336\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[0.1,1,3;1,0.2,3;5,2,0.3]\n",
+ "b=[12;2;4]\n",
+ "A\\b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ba =\n",
+ "\n",
+ " 120.000\n",
+ " 10.000\n",
+ " 13.333\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0 1 3\n",
+ " 1 0 3\n",
+ " 5 2 0\n",
+ "\n",
+ "sA =\n",
+ "\n",
+ " 0.00000 10.00000 30.00000\n",
+ " 5.00000 0.00000 15.00000\n",
+ " 16.66667 6.66667 0.00000\n",
+ "\n",
+ "x1 =\n",
+ "\n",
+ " 120.000\n",
+ " 10.000\n",
+ " 13.333\n",
+ "\n",
+ "x2 =\n",
+ "\n",
+ " -380.00\n",
+ " -790.00\n",
+ " -2053.33\n",
+ "\n",
+ "x3 =\n",
+ "\n",
+ " 6.9620e+04\n",
+ " 3.2710e+04\n",
+ " 1.1613e+04\n",
+ "\n",
+ "solution is not converging to [-2.93,9.19,1.03]\n"
+ ]
+ }
+ ],
+ "source": [
+ "ba=b./diag(A) % or ba=b./[A(1,1);A(2,2);A(3,3)]\n",
+ "sA=A-diag(diag(A)) % A with zeros on diagonal\n",
+ "sA(1,:)=sA(1,:)/A(1,1);\n",
+ "sA(2,:)=sA(2,:)/A(2,2);\n",
+ "sA(3,:)=sA(3,:)/A(3,3)\n",
+ "x0=[0;0;0];\n",
+ "x1=ba-sA*x0\n",
+ "x2=ba-sA*x1\n",
+ "x3=ba-sA*x2\n",
+ "fprintf('solution is not converging to [-2.93,9.19,1.03]\\n')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Gauss-Seidel with Relaxation\n",
+ "\n",
+ "In order to force the solution to converge faster, we can introduce a relaxation term $\\lambda$. \n",
+ "\n",
+ "where the new x values are weighted between the old and new:\n",
+ "\n",
+ "$x^{i}=\\lambda x^{i}+(1-\\lambda)x^{i-1}$\n",
+ "\n",
+ "after solving for x, lambda weights the current approximation with the previous approximation for the updated x\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "A =\n",
+ "\n",
+ " 3.00000 -0.10000 -0.20000\n",
+ " 0.10000 7.00000 -0.30000\n",
+ " 0.30000 -0.20000 10.00000\n",
+ "\n",
+ "b =\n",
+ "\n",
+ " 7.8500\n",
+ " -19.3000\n",
+ " 71.4000\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% rearrange A and b\n",
+ "A=[3 -0.1 -0.2;0.1 7 -0.3;0.3 -0.2 10]\n",
+ "b=[7.85;-19.3;71.4]\n",
+ "\n",
+ "iters=zeros(100,1);\n",
+ "for i=1:100\n",
+ " lambda=2/100*i;\n",
+ " [x,ea,iters(i)]=Jacobi_rel(A,b,lambda);\n",
+ "end\n",
+ "plot([1:100]*2/100,iters) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 107,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "l = 0.99158\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "l=fminbnd(@(l) lambda_fcn(A,b,l),0.5,1.5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 108,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A\\b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 109,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 1.8289e-07\n",
+ " 2.1984e-08\n",
+ " 2.3864e-08\n",
+ "\n",
+ "iter = 8\n",
+ "x =\n",
+ "\n",
+ " 3.0000\n",
+ " -2.5000\n",
+ " 7.0000\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 1.9130e-08\n",
+ " 7.6449e-08\n",
+ " 3.3378e-08\n",
+ "\n",
+ "iter = 8\n"
+ ]
+ }
+ ],
+ "source": [
+ "[x,ea,iter]=Jacobi_rel(A,b,l,0.000001)\n",
+ "[x,ea,iter]=Jacobi_rel(A,b,1,0.000001)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Nonlinear Systems\n",
+ "\n",
+ "Consider two simultaneous nonlinear equations with two unknowns:\n",
+ "\n",
+ "$x_{1}^{2}+x_{1}x_{2}=10$\n",
+ "\n",
+ "$x_{2}+3x_{1}x_{2}^{2}=57$\n",
+ "\n",
+ "Graphically, we are looking for the solution:\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x11=linspace(0.5,3);\n",
+ "x12=(10-x11.^2)./x11;\n",
+ "\n",
+ "x22=linspace(2,8);\n",
+ "x21=(57-x22).*x22.^-2/3;\n",
+ "\n",
+ "plot(x11,x12,x21,x22)\n",
+ "% Solution at x_1=2, x_2=3\n",
+ "hold on;\n",
+ "plot(2,3,'o')\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Newton-Raphson part II\n",
+ "\n",
+ "Remember the first order approximation for the next point in a function is:\n",
+ "\n",
+ "$f(x_{i+1})=f(x_{i})+(x_{i+1}-x_{i})f'(x_{i})$\n",
+ "\n",
+ "then, $f(x_{i+1})=0$ so we are left with:\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n",
+ "\n",
+ "We can use the same formula, but now we have multiple dimensions so we need to determine the Jacobian\n",
+ "\n",
+ "$[J]=\\left[ \\begin{array}{cccc}\n",
+ "\\frac{\\partial f_{1,i}}{\\partial x_{1}} & \\frac{\\partial f_{1,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{1,i}}{\\partial x_{n}} \\\\\n",
+ "\\frac{\\partial f_{2,i}}{\\partial x_{1}} & \\frac{\\partial f_{2,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{2,i}}{\\partial x_{n}} \\\\\n",
+ "\\vdots & \\vdots & & \\vdots \\\\\n",
+ "\\frac{\\partial f_{n,i}}{\\partial x_{1}} & \\frac{\\partial f_{n,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{n,i}}{\\partial x_{n}} \\\\\n",
+ "\\end{array} \\right]$\n",
+ "\n",
+ "$\\left[ \\begin{array}{c}\n",
+ "f_{1,i+1} \\\\\n",
+ "f_{2,i+1} \\\\\n",
+ "\\vdots \\\\\n",
+ "f_{n,i+1}\\end{array} \\right]=\n",
+ "\\left[ \\begin{array}{c}\n",
+ "f_{1,i} \\\\\n",
+ "f_{2,i} \\\\\n",
+ "\\vdots \\\\\n",
+ "f_{n,i}\\end{array} \\right]+\n",
+ "\\left[ \\begin{array}{cccc}\n",
+ "\\frac{\\partial f_{1,i}}{\\partial x_{1}} & \\frac{\\partial f_{1,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{1,i}}{\\partial x_{n}} \\\\\n",
+ "\\frac{\\partial f_{2,i}}{\\partial x_{1}} & \\frac{\\partial f_{2,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{2,i}}{\\partial x_{n}} \\\\\n",
+ "\\vdots & \\vdots & & \\vdots \\\\\n",
+ "\\frac{\\partial f_{n,i}}{\\partial x_{1}} & \\frac{\\partial f_{n,i}}{\\partial x_{2}} & \n",
+ "\\cdots & \\frac{\\partial f_{n,i}}{\\partial x_{n}} \\\\\n",
+ "\\end{array} \\right]\n",
+ "\\left( \\left[ \\begin{array}{c}\n",
+ "x_{i+1} \\\\\n",
+ "x_{i+1} \\\\\n",
+ "\\vdots \\\\\n",
+ "x_{i+1}\\end{array} \\right]-\n",
+ "\\left[ \\begin{array}{c}\n",
+ "x_{1,i} \\\\\n",
+ "x_{2,i} \\\\\n",
+ "\\vdots \\\\\n",
+ "x_{n,i}\\end{array} \\right]\\right)$\n",
+ "\n",
+ "### Solution is again in the form Ax=b\n",
+ "\n",
+ "$[J]([x_{i+1}]-[x_{i}])=-[f]$\n",
+ "\n",
+ "so\n",
+ "\n",
+ "$[x_{i+1}]= [x_{i}]-[J]^{-1}[f]$\n",
+ "\n",
+ "## Example of Jacobian calculation\n",
+ "\n",
+ "### Nonlinear springs supporting two masses in series\n",
+ "\n",
+ "Two springs are connected to two masses, with $m_1$=1 kg and $m_{2}$=2 kg. The springs are identical, but they have nonlinear spring constants, of $k_1$=100 N/m and $k_2$=-10 N/m\n",
+ "\n",
+ "We want to solve for the final position of the masses ($x_1$ and $x_2$)\n",
+ "\n",
+ "$m_{1}g+k_{1}(x_{2}-x_{1})+k_{2}(x_{2}-x_{1})^{2}+k_{1}x_{1}+k_{2}x_{1}^{2}=0$\n",
+ "\n",
+ "$m_{2}g-k_{1}(x_{2}-x_{1})-k_{2}(x_2-x_1)^{2}=0$\n",
+ "\n",
+ "$J(1,1)=\\frac{\\partial f_{1}}{\\partial x_{1}}=-k_{1}-2k_{2}(x_{2}-x_{1})+k_{1}+2k_{2}x_{1}$\n",
+ "\n",
+ "$J(1,2)=\\frac{\\partial f_1}{\\partial x_{2}}=k_{1}+2k_{2}(x_{2}-x_{1})$\n",
+ "\n",
+ "$J(2,1)=\\frac{\\partial f_2}{\\partial x_{1}}=k_{1}+2k_{2}(x_{2}-x_{1})$\n",
+ "\n",
+ "$J(2,2)=\\frac{\\partial f_2}{\\partial x_{2}}=-k_{1}-2k_{2}(x_{2}-x_{1})$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "m1=1; % kg \n",
+ "m2=2; % kg\n",
+ "k1=100; % N/m\n",
+ "k2=-10; % N/m^2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "function [f,J]=mass_spring(x)\n",
+ " % Function to calculate function values f1 and f2 as well as Jacobian \n",
+ " % for 2 masses and 2 identical nonlinear springs\n",
+ " m1=1; % kg \n",
+ " m2=2; % kg\n",
+ " k1=100; % N/m\n",
+ " k2=-10; % N/m^2\n",
+ " g=9.81; % m/s^2\n",
+ " x1=x(1);\n",
+ " x2=x(2);\n",
+ " J=[-k1-2*k2*(x2-x1)-k1-2*k2*x1,k1+2*k2*(x2-x1);\n",
+ " k1+2*k2*(x2-x1),-k1-2*k2*(x2-x1)];\n",
+ " f=[m1*g+k1*(x2-x1)+k2*(x2-x1).^2-k1*x1-k2*x1^2;\n",
+ " m2*g-k1*(x2-x1)-k2*(x2-x1).^2];\n",
+ "end\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "f =\n",
+ "\n",
+ " -190.19\n",
+ " 129.62\n",
+ "\n",
+ "J =\n",
+ "\n",
+ " -200 120\n",
+ " 120 -120\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "[f,J]=mass_spring([1,0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x1 =\n",
+ "\n",
+ " -1.5142\n",
+ " -1.4341\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 2.9812\n",
+ " 2.3946\n",
+ "\n",
+ "x2 =\n",
+ "\n",
+ " 0.049894\n",
+ " 0.248638\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 31.3492\n",
+ " 6.7678\n",
+ "\n",
+ "x3 =\n",
+ "\n",
+ " 0.29701\n",
+ " 0.49722\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 0.83201\n",
+ " 0.49995\n",
+ "\n",
+ "x =\n",
+ "\n",
+ " 0.29701\n",
+ " 0.49722\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 0.021392\n",
+ " 0.012890\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 1.4786e-05\n",
+ " 8.9091e-06\n",
+ "\n",
+ "ea =\n",
+ "\n",
+ " 7.0642e-12\n",
+ " 4.2565e-12\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x0=[3;2];\n",
+ "[f0,J0]=mass_spring(x0);\n",
+ "x1=x0-J0\\f0\n",
+ "ea=(x1-x0)./x1\n",
+ "[f1,J1]=mass_spring(x1);\n",
+ "x2=x1-J1\\f1\n",
+ "ea=(x2-x1)./x2\n",
+ "[f2,J2]=mass_spring(x2);\n",
+ "x3=x2-J2\\f2\n",
+ "ea=(x3-x2)./x3\n",
+ "x=x3\n",
+ "for i=1:3\n",
+ " xold=x;\n",
+ " [f,J]=mass_spring(x);\n",
+ " x=x-J\\f;\n",
+ " ea=(x-xold)./x\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x =\n",
+ "\n",
+ " 0.30351\n",
+ " 0.50372\n",
+ "\n",
+ "X0 =\n",
+ "\n",
+ " 0.30351\n",
+ " 0.50372\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "x\n",
+ "X0=fsolve(@(x) mass_spring(x),[3;5])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[X,Y]=meshgrid(linspace(0,10,20),linspace(0,10,20));\n",
+ "[N,M]=size(X);\n",
+ "F=zeros(size(X));\n",
+ "for i=1:N\n",
+ " for j=1:M\n",
+ " [f,~]=mass_spring([X(i,j),Y(i,j)]);\n",
+ " F1(i,j)=f(1);\n",
+ " F2(i,j)=f(2);\n",
+ " end\n",
+ "end\n",
+ "mesh(X,Y,F1)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')\n",
+ "colorbar()\n",
+ "figure()\n",
+ "mesh(X,Y,F2)\n",
+ "xlabel('x_1')\n",
+ "ylabel('x_2')\n",
+ "colorbar()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/12_matrix_condition/GS_rel.m b/12_matrix_condition/GS_rel.m
new file mode 100644
index 0000000..4a6daf4
--- /dev/null
+++ b/12_matrix_condition/GS_rel.m
@@ -0,0 +1,36 @@
+function [x,ea,iter] = GS_rel(A,b,lambda,es,maxit)
+% GaussSeidel: Gauss Seidel method
+% x = GaussSeidel(A,b): Gauss Seidel without relaxation
+% input:
+% A = coefficient matrix
+% b = right hand side vector
+% es = stop criterion (default = 0.00001%)
+% maxit = max iterations (default = 50)
+% output:
+% x = solution vector
+if nargin<3,error('at least 2 input arguments required'),end
+if nargin<5|isempty(maxit),maxit=50;end
+if nargin<4|isempty(es),es=0.00001;end
+[m,n] = size(A);
+if m~=n, error('Matrix A must be square'); end
+C = A-diag(diag(A));
+x=zeros(n,1);
+for i = 1:n
+ C(i,1:n) = C(i,1:n)/A(i,i);
+end
+
+d = b./diag(A);
+
+iter = 0;
+while (1)
+ xold = x;
+ for i = 1:n
+ x(i) = d(i)-C(i,:)*x;
+ x(i) = lambda*x(i)+(1-lambda)*xold(i);
+ if x(i) ~= 0
+ ea(i) = abs((x(i) - xold(i))/x(i)) * 100;
+ end
+ end
+ iter = iter+1;
+ if max(ea)<=es | iter >= maxit, break, end
+end
diff --git a/12_matrix_condition/GaussSeidel.m b/12_matrix_condition/GaussSeidel.m
new file mode 100644
index 0000000..2be52e1
--- /dev/null
+++ b/12_matrix_condition/GaussSeidel.m
@@ -0,0 +1,35 @@
+function x = GaussSeidel(A,b,es,maxit)
+% GaussSeidel: Gauss Seidel method
+% x = GaussSeidel(A,b): Gauss Seidel without relaxation
+% input:
+% A = coefficient matrix
+% b = right hand side vector
+% es = stop criterion (default = 0.00001%)
+% maxit = max iterations (default = 50)
+% output:
+% x = solution vector
+if nargin<2,error('at least 2 input arguments required'),end
+if nargin<4|isempty(maxit),maxit=50;end
+if nargin<3|isempty(es),es=0.00001;end
+[m,n] = size(A);
+if m~=n, error('Matrix A must be square'); end
+C = A-diag(diag(A));
+x=zeros(n,1);
+for i = 1:n
+ C(i,1:n) = C(i,1:n)/A(i,i);
+end
+
+d = b./diag(A);
+
+iter = 0;
+while (1)
+ xold = x;
+ for i = 1:n
+ x(i) = d(i)-C(i,:)*x;
+ if x(i) ~= 0
+ ea(i) = abs((x(i) - xold(i))/x(i)) * 100;
+ end
+ end
+ iter = iter+1;
+ if max(ea)<=es | iter >= maxit, break, end
+end
diff --git a/12_matrix_condition/Jacobi.m b/12_matrix_condition/Jacobi.m
new file mode 100644
index 0000000..8a7b4ae
--- /dev/null
+++ b/12_matrix_condition/Jacobi.m
@@ -0,0 +1,39 @@
+function x = Jacobi(A,b,es,maxit)
+% GaussSeidel: Gauss Seidel method
+% x = GaussSeidel(A,b): Gauss Seidel without relaxation
+% input:
+% A = coefficient matrix
+% b = right hand side vector
+% es = stop criterion (default = 0.00001%)
+% maxit = max iterations (default = 50)
+% output:
+% x = solution vector
+if nargin<2,error('at least 2 input arguments required'),end
+if nargin<4|isempty(maxit),maxit=50;end
+if nargin<3|isempty(es),es=0.00001;end
+[m,n] = size(A);
+if m~=n, error('Matrix A must be square'); end
+C = A-diag(diag(A));
+x=zeros(n,1);
+for i = 1:n
+ C(i,1:n) = C(i,1:n)/A(i,i);
+end
+
+d = b./diag(A);
+
+iter = 0;
+while (1)
+ xold = x;
+ x = d-C*x;
+ % if any values of x are zero, we add 1 to denominator so error is well-behaved
+ i_zero=find(x==0);
+ i=find(x~=0);
+ if length(i_zero)>0
+ ea(i_zero)=abs((x-xold)./(1+x)*100);
+ ea(i) = abs((x(i) - xold(i))./x(i)) * 100;
+ else
+ ea = abs((x - xold)./x) * 100;
+ end
+ iter = iter+1;
+ if max(ea)<=es | iter >= maxit, break, end
+end
diff --git a/12_matrix_condition/Jacobi_rel.m b/12_matrix_condition/Jacobi_rel.m
new file mode 100644
index 0000000..5cdec33
--- /dev/null
+++ b/12_matrix_condition/Jacobi_rel.m
@@ -0,0 +1,41 @@
+function [x,ea,iter]= Jacobi_rel(A,b,lambda,es,maxit)
+% GaussSeidel: Gauss Seidel method
+% x = GaussSeidel(A,b): Gauss Seidel without relaxation
+% input:
+% A = coefficient matrix
+% b = right hand side vector
+% es = stop criterion (default = 0.00001%)
+% maxit = max iterations (default = 50)
+% output:
+% x = solution vector
+if nargin<3,error('at least 2 input arguments required'),end
+if nargin<5|isempty(maxit),maxit=50;end
+if nargin<4|isempty(es),es=0.00001;end
+[m,n] = size(A);
+if m~=n, error('Matrix A must be square'); end
+C = A-diag(diag(A));
+x=zeros(n,1);
+for i = 1:n
+ C(i,1:n) = C(i,1:n)/A(i,i);
+end
+
+d = b./diag(A);
+
+iter = 0;
+while (1)
+ xold = x;
+ x = d-C*x;
+ % Add relaxation parameter lambda to current iteration
+ x = lambda*x+(1-lambda)*xold;
+ % if any values of x are zero, we add 1 to denominator so error is well-behaved
+ i_zero=find(x==0);
+ i=find(x~=0);
+ if length(i_zero)>0
+ ea(i_zero)=abs((x-xold)./(1+x)*100);
+ ea(i) = abs((x(i) - xold(i))./x(i)) * 100;
+ else
+ ea = abs((x - xold)./x) * 100;
+ end
+ iter = iter+1;
+ if max(ea)<=es | iter >= maxit, break, end
+end
diff --git a/12_matrix_condition/LU_naive.m b/12_matrix_condition/LU_naive.m
new file mode 100644
index 0000000..92efde6
--- /dev/null
+++ b/12_matrix_condition/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/12_matrix_condition/efficient_soln.png b/12_matrix_condition/efficient_soln.png
new file mode 100644
index 0000000..ef24ece
Binary files /dev/null and b/12_matrix_condition/efficient_soln.png differ
diff --git a/12_matrix_condition/gp_image_01.png b/12_matrix_condition/gp_image_01.png
new file mode 100644
index 0000000..ef291b5
Binary files /dev/null and b/12_matrix_condition/gp_image_01.png differ
diff --git a/12_matrix_condition/lambda_fcn.m b/12_matrix_condition/lambda_fcn.m
new file mode 100644
index 0000000..435f9e9
--- /dev/null
+++ b/12_matrix_condition/lambda_fcn.m
@@ -0,0 +1,8 @@
+function iters = lambda_fcn(A,b,lambda)
+ % function to minimize the number of iterations for a given Ax=b solution
+ % using default Jacobi_rel parameters of es=0.00001 and maxit=50
+
+ [x,ea,iters]= Jacobi_rel(A,b,lambda,1e-8);
+end
+
+
diff --git a/12_matrix_condition/norm_A.png b/12_matrix_condition/norm_A.png
new file mode 100644
index 0000000..5f2d273
Binary files /dev/null and b/12_matrix_condition/norm_A.png differ
diff --git a/12_matrix_condition/octave-workspace b/12_matrix_condition/octave-workspace
new file mode 100644
index 0000000..8c437bb
Binary files /dev/null and b/12_matrix_condition/octave-workspace differ
diff --git a/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb b/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb
new file mode 100644
index 0000000..f6b2b95
--- /dev/null
+++ b/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb
@@ -0,0 +1,590 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Eigenvalues\n",
+ "\n",
+ "Eigenvalues and eigen vectors are the solution to the set of equations where\n",
+ "\n",
+ "$Ax=\\lambda x$\n",
+ "\n",
+ "or \n",
+ "\n",
+ "$A-I \\lambda=0$\n",
+ "\n",
+ "Where A is the description of the system and I is the identity matrix with the same dimensions as A and $\\lambda$ is an eigenvalue of A. \n",
+ "\n",
+ "These problems are seen in a number of engineering practices:\n",
+ "\n",
+ "1. Determining vibrational modes in structural devices\n",
+ "\n",
+ "2. Material science - vibrational modes of crystal lattices (phonons)\n",
+ "\n",
+ "3. [Google searches - http://www.rose-hulman.edu/~bryan/googleFinalVersionFixed.pdf](http://www.rose-hulman.edu/~bryan/googleFinalVersionFixed.pdf)\n",
+ "\n",
+ "4. Quantum mechanics - many solutions are eigenvalue problems\n",
+ "\n",
+ "5. Solid mechanics, principle stresses and principle stress directions are eigenvalues and eigenvectors\n",
+ "\n",
+ "One way of determining the eigenvalues is taking the determinant:\n",
+ "\n",
+ "$|A-\\lambda I|=0$\n",
+ "\n",
+ "This will result in an $n^{th}$-order polynomial where A is $n \\times n$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Take, A\n",
+ "\n",
+ "$A=\\left[\\begin{array}{ccc}\n",
+ "2 & 1 & 0 \\\\\n",
+ "1 & 3 & 1 \\\\\n",
+ "0 & 1 & 4 \\end{array}\\right]$\n",
+ "\n",
+ "$|A-\\lambda I| = \\left|\\begin{array}{ccc}\n",
+ "2-\\lambda & 1 & 0 \\\\\n",
+ "1 & 3-\\lambda & 1 \\\\\n",
+ "0 & 1 & 4-\\lambda \\end{array}\\right|=0$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$(2-\\lambda)(3-\\lambda)(4-\\lambda)-(4-\\lambda)-(2-\\lambda)=0$\n",
+ "\n",
+ "$-\\lambda^{3}+9\\lambda^{2}-24\\lambda+18 =0$\n",
+ "\n",
+ "$\\lambda = 3,~\\sqrt{3}+3,-\\sqrt{3}+3$\n",
+ "\n",
+ "in Matlab/Octave:\n",
+ "\n",
+ "```matlab\n",
+ "A=[2,1,0;1,3,1;0,1,4];\n",
+ "pA=poly(A);\n",
+ "lambda = roots(pA)\n",
+ "```\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "lambda =\n",
+ "\n",
+ " 4.7321\n",
+ " 3.0000\n",
+ " 1.2679\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[2,1,0;1,3,1;0,1,4];\n",
+ "pA=poly(A); % characteristic polynomial of A, e.g. l^3-9*l^2+24*l-18=0\n",
+ "lambda = roots(pA)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Applications of Eigenvalue analysis\n",
+ "\n",
+ "Consider the 2-mass, 3-spring system shown below\n",
+ "\n",
+ "\n",
+ "\n",
+ "It might not be immediately obvious, but there are two resonant frequencies for these masses connected in series. \n",
+ "\n",
+ "Take the two FBD solutions:\n",
+ "\n",
+ "$m_{1}\\frac{d^{2}x_{1}}{dt^{2}}=-kx_{1}+k(x_{2}-x_{1})$\n",
+ "\n",
+ "$m_{2}\\frac{d^{2}x_{2}}{dt^{2}}=-k(x_{2}-x_{1})-kx_{2}$\n",
+ "\n",
+ "we know that $x_{i}(t)\\propto sin(\\omega t)$ so we can substitute\n",
+ "\n",
+ "$x_{i}=X_{i}sin(\\omega t)$\n",
+ "\n",
+ "$-m_{1}X_{1}\\omega^{2}sin(\\omega t) = -kX_{1}sin(\\omega t) +k(X_{2}-X_{1})sin(\\omega t)$\n",
+ "\n",
+ "$-m_{2}X_{2}\\omega^{2}sin(\\omega t) = -kX_{2}sin(\\omega t) - k(X_{2}-X_{1})sin(\\omega t)$\n",
+ "\n",
+ "where $X_{1}$ and $X_{2}$ are the amplitude of oscillations and $\\omega$ is the frequency of oscillations. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "now, \n",
+ "\n",
+ "$\\left(\\frac{2k}{m_{1}}-\\omega^{2}\\right)X_{1}-\\frac{k}{m_{1}}X_{2} = 0$\n",
+ "\n",
+ "\n",
+ "$-\\frac{k}{m_{2}}X_{1}+\\left(\\frac{2k}{m_{2}}-\\omega^{2}\\right)X_{2} = 0$\n",
+ "\n",
+ "or\n",
+ "\n",
+ "$|K-\\lambda I| = \\left|\\begin{array}{cc}\n",
+ "\\left(\\frac{2k}{m_{1}}-\\omega^{2}\\right) & -\\frac{k}{m_{1}} \\\\\n",
+ "-\\frac{k}{m_{2}} & \\left(\\frac{2k}{m_{2}}-\\omega^{2}\\right)\n",
+ "\\end{array}\\right|=0$\n",
+ "\n",
+ "where $\\lambda = \\omega^{2}$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'eig' is a built-in function from the file libinterp/corefcn/eig.cc\n",
+ "\n",
+ " -- Built-in Function: LAMBDA = eig (A)\n",
+ " -- Built-in Function: LAMBDA = eig (A, B)\n",
+ " -- Built-in Function: [V, LAMBDA] = eig (A)\n",
+ " -- Built-in Function: [V, LAMBDA] = eig (A, B)\n",
+ " Compute the eigenvalues (and optionally the eigenvectors) of a\n",
+ " matrix or a pair of matrices\n",
+ "\n",
+ " The algorithm used depends on whether there are one or two input\n",
+ " matrices, if they are real or complex, and if they are symmetric\n",
+ " (Hermitian if complex) or non-symmetric.\n",
+ "\n",
+ " The eigenvalues returned by 'eig' are not ordered.\n",
+ "\n",
+ " See also: eigs, 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 eig"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 10 -5\n",
+ " -5 10\n",
+ "\n",
+ "v =\n",
+ "\n",
+ " -0.70711 -0.70711\n",
+ " -0.70711 0.70711\n",
+ "\n",
+ "e =\n",
+ "\n",
+ "Diagonal Matrix\n",
+ "\n",
+ " 5 0\n",
+ " 0 15\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "m=40; % mass in kg\n",
+ "k=200; % spring constant in N/m\n",
+ "\n",
+ "K=[2*k/m,-k/m;-k/m,2*k/m]\n",
+ "\n",
+ "[v,e]=eig(K)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " -10.607\n",
+ " 10.607\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -10.607\n",
+ " 10.607\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K*v(:,2)\n",
+ "e(2,2)*v(:,2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### $m_{1}=m_{2}=$40 kg and $k_{1}=k_{2}=k_{3}=$200 N/m\n",
+ "\n",
+ "\n",
+ "### $m_{1}=$40 kg, $m_{2}=$50 kg, $k_{1}=$200 N/m, and $k_{2}=k_{3}=$100 N/m\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solid Mechanics\n",
+ "\n",
+ "Many times you want to know the \"principle\" components of stress or strain. \n",
+ "\n",
+ "stress and strain are second order tensors:\n",
+ "\n",
+ "$\\sigma_{ij} =\\frac{F_{i}}{A_{j}}$ (engineering stress)\n",
+ "\n",
+ "$\\epsilon_{ij}=\\frac{\\Delta l_{i}}{l_{j}}$ (engineering strain)\n",
+ "\n",
+ "Imagine you can apply a force in two directions normal to each other, could you create a shear stress? \n",
+ "\n",
+ "\n",
+ "\n",
+ "Let's try to create the stress state of 10 MPa shear stress with two normal stresses. \n",
+ "\n",
+ "that means, $\\sigma_{12}=\\sigma_{21}=$10 MPa, and $\\sigma_{11}=\\sigma_{22}=\\sigma_{33}=\\sigma_{23}=\\sigma_{13}=0$ MPa\n",
+ "\n",
+ "in matrix form:\n",
+ "\n",
+ "$[\\sigma]=\\left[\\begin{array}{ccc}\n",
+ "0 & 10 & 0 \\\\\n",
+ "10 & 0 & 0 \\\\\n",
+ "0 & 0 & 0 \\end{array} \\right]$ MPa"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "v =\n",
+ "\n",
+ " -0.70711 0.00000 0.70711\n",
+ " 0.70711 0.00000 0.70711\n",
+ " 0.00000 1.00000 0.00000\n",
+ "\n",
+ "e =\n",
+ "\n",
+ "Diagonal Matrix\n",
+ "\n",
+ " -10 0 0\n",
+ " 0 0 0\n",
+ " 0 0 10\n",
+ "\n",
+ "v1 =\n",
+ "\n",
+ " 7.07107\n",
+ " -7.07107\n",
+ " 0.00000\n",
+ "\n",
+ "v2 =\n",
+ "\n",
+ " 7.07107\n",
+ " 7.07107\n",
+ " 0.00000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "s=[0,10,0;10,0,0;0,0,0];\n",
+ "[v,e]=eig(s)\n",
+ "v1=s*v(:,1)\n",
+ "v2=s*v(:,3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "v1=s*v(:,1);\n",
+ "v2=s*v(:,3);\n",
+ "%quiver([0,0],[0,0],v(1,[1,3]),v(2,[1,3]))\n",
+ "quiver([0,0],[0,0],[v1(1),v2(1)],[v1(2),v2(2)])\n",
+ "axis([-10,10,-10,10])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "e =\n",
+ "\n",
+ " -10\n",
+ " 0\n",
+ " 10\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=eig(s)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/13_eigenvalues/13_eigenvalues.ipynb b/13_eigenvalues/13_eigenvalues.ipynb
new file mode 100644
index 0000000..f6b2b95
--- /dev/null
+++ b/13_eigenvalues/13_eigenvalues.ipynb
@@ -0,0 +1,590 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Eigenvalues\n",
+ "\n",
+ "Eigenvalues and eigen vectors are the solution to the set of equations where\n",
+ "\n",
+ "$Ax=\\lambda x$\n",
+ "\n",
+ "or \n",
+ "\n",
+ "$A-I \\lambda=0$\n",
+ "\n",
+ "Where A is the description of the system and I is the identity matrix with the same dimensions as A and $\\lambda$ is an eigenvalue of A. \n",
+ "\n",
+ "These problems are seen in a number of engineering practices:\n",
+ "\n",
+ "1. Determining vibrational modes in structural devices\n",
+ "\n",
+ "2. Material science - vibrational modes of crystal lattices (phonons)\n",
+ "\n",
+ "3. [Google searches - http://www.rose-hulman.edu/~bryan/googleFinalVersionFixed.pdf](http://www.rose-hulman.edu/~bryan/googleFinalVersionFixed.pdf)\n",
+ "\n",
+ "4. Quantum mechanics - many solutions are eigenvalue problems\n",
+ "\n",
+ "5. Solid mechanics, principle stresses and principle stress directions are eigenvalues and eigenvectors\n",
+ "\n",
+ "One way of determining the eigenvalues is taking the determinant:\n",
+ "\n",
+ "$|A-\\lambda I|=0$\n",
+ "\n",
+ "This will result in an $n^{th}$-order polynomial where A is $n \\times n$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Take, A\n",
+ "\n",
+ "$A=\\left[\\begin{array}{ccc}\n",
+ "2 & 1 & 0 \\\\\n",
+ "1 & 3 & 1 \\\\\n",
+ "0 & 1 & 4 \\end{array}\\right]$\n",
+ "\n",
+ "$|A-\\lambda I| = \\left|\\begin{array}{ccc}\n",
+ "2-\\lambda & 1 & 0 \\\\\n",
+ "1 & 3-\\lambda & 1 \\\\\n",
+ "0 & 1 & 4-\\lambda \\end{array}\\right|=0$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$(2-\\lambda)(3-\\lambda)(4-\\lambda)-(4-\\lambda)-(2-\\lambda)=0$\n",
+ "\n",
+ "$-\\lambda^{3}+9\\lambda^{2}-24\\lambda+18 =0$\n",
+ "\n",
+ "$\\lambda = 3,~\\sqrt{3}+3,-\\sqrt{3}+3$\n",
+ "\n",
+ "in Matlab/Octave:\n",
+ "\n",
+ "```matlab\n",
+ "A=[2,1,0;1,3,1;0,1,4];\n",
+ "pA=poly(A);\n",
+ "lambda = roots(pA)\n",
+ "```\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "lambda =\n",
+ "\n",
+ " 4.7321\n",
+ " 3.0000\n",
+ " 1.2679\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=[2,1,0;1,3,1;0,1,4];\n",
+ "pA=poly(A); % characteristic polynomial of A, e.g. l^3-9*l^2+24*l-18=0\n",
+ "lambda = roots(pA)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Applications of Eigenvalue analysis\n",
+ "\n",
+ "Consider the 2-mass, 3-spring system shown below\n",
+ "\n",
+ "\n",
+ "\n",
+ "It might not be immediately obvious, but there are two resonant frequencies for these masses connected in series. \n",
+ "\n",
+ "Take the two FBD solutions:\n",
+ "\n",
+ "$m_{1}\\frac{d^{2}x_{1}}{dt^{2}}=-kx_{1}+k(x_{2}-x_{1})$\n",
+ "\n",
+ "$m_{2}\\frac{d^{2}x_{2}}{dt^{2}}=-k(x_{2}-x_{1})-kx_{2}$\n",
+ "\n",
+ "we know that $x_{i}(t)\\propto sin(\\omega t)$ so we can substitute\n",
+ "\n",
+ "$x_{i}=X_{i}sin(\\omega t)$\n",
+ "\n",
+ "$-m_{1}X_{1}\\omega^{2}sin(\\omega t) = -kX_{1}sin(\\omega t) +k(X_{2}-X_{1})sin(\\omega t)$\n",
+ "\n",
+ "$-m_{2}X_{2}\\omega^{2}sin(\\omega t) = -kX_{2}sin(\\omega t) - k(X_{2}-X_{1})sin(\\omega t)$\n",
+ "\n",
+ "where $X_{1}$ and $X_{2}$ are the amplitude of oscillations and $\\omega$ is the frequency of oscillations. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "now, \n",
+ "\n",
+ "$\\left(\\frac{2k}{m_{1}}-\\omega^{2}\\right)X_{1}-\\frac{k}{m_{1}}X_{2} = 0$\n",
+ "\n",
+ "\n",
+ "$-\\frac{k}{m_{2}}X_{1}+\\left(\\frac{2k}{m_{2}}-\\omega^{2}\\right)X_{2} = 0$\n",
+ "\n",
+ "or\n",
+ "\n",
+ "$|K-\\lambda I| = \\left|\\begin{array}{cc}\n",
+ "\\left(\\frac{2k}{m_{1}}-\\omega^{2}\\right) & -\\frac{k}{m_{1}} \\\\\n",
+ "-\\frac{k}{m_{2}} & \\left(\\frac{2k}{m_{2}}-\\omega^{2}\\right)\n",
+ "\\end{array}\\right|=0$\n",
+ "\n",
+ "where $\\lambda = \\omega^{2}$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'eig' is a built-in function from the file libinterp/corefcn/eig.cc\n",
+ "\n",
+ " -- Built-in Function: LAMBDA = eig (A)\n",
+ " -- Built-in Function: LAMBDA = eig (A, B)\n",
+ " -- Built-in Function: [V, LAMBDA] = eig (A)\n",
+ " -- Built-in Function: [V, LAMBDA] = eig (A, B)\n",
+ " Compute the eigenvalues (and optionally the eigenvectors) of a\n",
+ " matrix or a pair of matrices\n",
+ "\n",
+ " The algorithm used depends on whether there are one or two input\n",
+ " matrices, if they are real or complex, and if they are symmetric\n",
+ " (Hermitian if complex) or non-symmetric.\n",
+ "\n",
+ " The eigenvalues returned by 'eig' are not ordered.\n",
+ "\n",
+ " See also: eigs, 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 eig"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "K =\n",
+ "\n",
+ " 10 -5\n",
+ " -5 10\n",
+ "\n",
+ "v =\n",
+ "\n",
+ " -0.70711 -0.70711\n",
+ " -0.70711 0.70711\n",
+ "\n",
+ "e =\n",
+ "\n",
+ "Diagonal Matrix\n",
+ "\n",
+ " 5 0\n",
+ " 0 15\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "m=40; % mass in kg\n",
+ "k=200; % spring constant in N/m\n",
+ "\n",
+ "K=[2*k/m,-k/m;-k/m,2*k/m]\n",
+ "\n",
+ "[v,e]=eig(K)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " -10.607\n",
+ " 10.607\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " -10.607\n",
+ " 10.607\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "K*v(:,2)\n",
+ "e(2,2)*v(:,2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### $m_{1}=m_{2}=$40 kg and $k_{1}=k_{2}=k_{3}=$200 N/m\n",
+ "\n",
+ "\n",
+ "### $m_{1}=$40 kg, $m_{2}=$50 kg, $k_{1}=$200 N/m, and $k_{2}=k_{3}=$100 N/m\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solid Mechanics\n",
+ "\n",
+ "Many times you want to know the \"principle\" components of stress or strain. \n",
+ "\n",
+ "stress and strain are second order tensors:\n",
+ "\n",
+ "$\\sigma_{ij} =\\frac{F_{i}}{A_{j}}$ (engineering stress)\n",
+ "\n",
+ "$\\epsilon_{ij}=\\frac{\\Delta l_{i}}{l_{j}}$ (engineering strain)\n",
+ "\n",
+ "Imagine you can apply a force in two directions normal to each other, could you create a shear stress? \n",
+ "\n",
+ "\n",
+ "\n",
+ "Let's try to create the stress state of 10 MPa shear stress with two normal stresses. \n",
+ "\n",
+ "that means, $\\sigma_{12}=\\sigma_{21}=$10 MPa, and $\\sigma_{11}=\\sigma_{22}=\\sigma_{33}=\\sigma_{23}=\\sigma_{13}=0$ MPa\n",
+ "\n",
+ "in matrix form:\n",
+ "\n",
+ "$[\\sigma]=\\left[\\begin{array}{ccc}\n",
+ "0 & 10 & 0 \\\\\n",
+ "10 & 0 & 0 \\\\\n",
+ "0 & 0 & 0 \\end{array} \\right]$ MPa"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "v =\n",
+ "\n",
+ " -0.70711 0.00000 0.70711\n",
+ " 0.70711 0.00000 0.70711\n",
+ " 0.00000 1.00000 0.00000\n",
+ "\n",
+ "e =\n",
+ "\n",
+ "Diagonal Matrix\n",
+ "\n",
+ " -10 0 0\n",
+ " 0 0 0\n",
+ " 0 0 10\n",
+ "\n",
+ "v1 =\n",
+ "\n",
+ " 7.07107\n",
+ " -7.07107\n",
+ " 0.00000\n",
+ "\n",
+ "v2 =\n",
+ "\n",
+ " 7.07107\n",
+ " 7.07107\n",
+ " 0.00000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "s=[0,10,0;10,0,0;0,0,0];\n",
+ "[v,e]=eig(s)\n",
+ "v1=s*v(:,1)\n",
+ "v2=s*v(:,3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "v1=s*v(:,1);\n",
+ "v2=s*v(:,3);\n",
+ "%quiver([0,0],[0,0],v(1,[1,3]),v(2,[1,3]))\n",
+ "quiver([0,0],[0,0],[v1(1),v2(1)],[v1(2),v2(2)])\n",
+ "axis([-10,10,-10,10])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "e =\n",
+ "\n",
+ " -10\n",
+ " 0\n",
+ " 10\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "e=eig(s)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Octave",
+ "language": "octave",
+ "name": "octave"
+ },
+ "language_info": {
+ "file_extension": ".m",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "octave",
+ "version": "0.19.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/13_eigenvalues/animate_eig.m b/13_eigenvalues/animate_eig.m
new file mode 100644
index 0000000..a4721a5
--- /dev/null
+++ b/13_eigenvalues/animate_eig.m
@@ -0,0 +1,56 @@
+function [v,e]=animate_eig(k,m)
+ % create a series of png's for use in an animation based on
+ % spring constants [k1,k2,k3]=k
+ % and
+ % masses
+ % [m1,m2]=m
+
+ % check inputs for m and k
+ if length(m)==1
+ m1=m;
+ m2=m;
+ else
+ m1=m(1);
+ m2=m(2);
+ end
+ if length(k)==1
+ k1=k; k2=k; k3=k;
+ else
+ k1=k(1);
+ k2=k(2);
+ k3=k(3);
+ end
+ setdefaults
+ K=[k1/m1+k2/m1,-k2/m1;-k3/m2,k2/m2+k3/m2];
+ [v,e]=eig(K);
+ w1=e(1,1); w2=e(2,2);
+ scale = 0.5; % the magnitude of oscillations is not important for vibrational modes
+ % just set scale to 1/2 for nice plots
+ % the eigenvector magnitude is independent of the solution
+ X11=v(1,1)*scale; X12=v(2,1)*scale;
+ X21=v(1,2)*scale; X22=v(2,2)*scale;
+ f11=@(t) X11*sin(w1*t); f12=@(t) X12*sin(w1*t);
+ f21=@(t) X21*sin(w2*t); f22=@(t) X22*sin(w2*t);
+ f=figure();
+ time=linspace(-1,4);
+ % create a loop to plot the position over time where mass 1 and 2 start at x=1 and x=2 m
+ % then for the next vibrational mode, mass 1 and 2 start at x=3 and x=4 m
+
+ for i=1:length(time)
+ t=time(i);
+ plot(f11(t)+1,0,'rs',f11(time+t)+1,-time,...
+ f12(t)+2,0,'rs',f12(time+t)+2,-time,...
+ f21(t)+3,0,'bo',f21(time+t)+3,-time,...
+ f22(t)+4,0,'bo',f22(time+t)+4,-time)
+ axis([0 6 -4 1])
+ title('Vibration Modes')
+ xlabel('position (m)')
+ ylabel('time (s)')
+ filename=sprintf('output/%05d.png',i);
+ % another option for saving a plot to a png is the 'print' command
+ print(filename)
+ end
+ % this is a system command to use ImageMagick's cli 'convert' to create an animated gif
+ % if you don't have ImageMagick installed, comment this next line
+ system ("convert -delay 10 -loop 0 output/*png eigenvalues.gif")
+end
diff --git a/13_eigenvalues/eig_200_100_40_20.gif b/13_eigenvalues/eig_200_100_40_20.gif
new file mode 100644
index 0000000..1b583a7
Binary files /dev/null and b/13_eigenvalues/eig_200_100_40_20.gif differ
diff --git a/13_eigenvalues/eig_200_100_40_50.gif b/13_eigenvalues/eig_200_100_40_50.gif
new file mode 100644
index 0000000..ba05f59
Binary files /dev/null and b/13_eigenvalues/eig_200_100_40_50.gif differ
diff --git a/13_eigenvalues/eig_200_40.gif b/13_eigenvalues/eig_200_40.gif
new file mode 100644
index 0000000..f29f27e
Binary files /dev/null and b/13_eigenvalues/eig_200_40.gif differ
diff --git a/13_eigenvalues/eigenvalues.gif b/13_eigenvalues/eigenvalues.gif
new file mode 100644
index 0000000..1286762
Binary files /dev/null and b/13_eigenvalues/eigenvalues.gif differ
diff --git a/13_eigenvalues/octave-workspace b/13_eigenvalues/octave-workspace
new file mode 100644
index 0000000..8c437bb
Binary files /dev/null and b/13_eigenvalues/octave-workspace differ
diff --git a/13_eigenvalues/output/00001.png b/13_eigenvalues/output/00001.png
new file mode 100644
index 0000000..44210b6
Binary files /dev/null and b/13_eigenvalues/output/00001.png differ
diff --git a/13_eigenvalues/output/00002.png b/13_eigenvalues/output/00002.png
new file mode 100644
index 0000000..31a03bb
Binary files /dev/null and b/13_eigenvalues/output/00002.png differ
diff --git a/13_eigenvalues/output/00003.png b/13_eigenvalues/output/00003.png
new file mode 100644
index 0000000..08cd507
Binary files /dev/null and b/13_eigenvalues/output/00003.png differ
diff --git a/13_eigenvalues/output/00004.png b/13_eigenvalues/output/00004.png
new file mode 100644
index 0000000..590a413
Binary files /dev/null and b/13_eigenvalues/output/00004.png differ
diff --git a/13_eigenvalues/output/00005.png b/13_eigenvalues/output/00005.png
new file mode 100644
index 0000000..240744a
Binary files /dev/null and b/13_eigenvalues/output/00005.png differ
diff --git a/13_eigenvalues/output/00006.png b/13_eigenvalues/output/00006.png
new file mode 100644
index 0000000..e31a3c7
Binary files /dev/null and b/13_eigenvalues/output/00006.png differ
diff --git a/13_eigenvalues/output/00007.png b/13_eigenvalues/output/00007.png
new file mode 100644
index 0000000..7ad26e6
Binary files /dev/null and b/13_eigenvalues/output/00007.png differ
diff --git a/13_eigenvalues/output/00008.png b/13_eigenvalues/output/00008.png
new file mode 100644
index 0000000..10bcf68
Binary files /dev/null and b/13_eigenvalues/output/00008.png differ
diff --git a/13_eigenvalues/output/00009.png b/13_eigenvalues/output/00009.png
new file mode 100644
index 0000000..53e6f1c
Binary files /dev/null and b/13_eigenvalues/output/00009.png differ
diff --git a/13_eigenvalues/output/00010.png b/13_eigenvalues/output/00010.png
new file mode 100644
index 0000000..c78c8ca
Binary files /dev/null and b/13_eigenvalues/output/00010.png differ
diff --git a/13_eigenvalues/output/00011.png b/13_eigenvalues/output/00011.png
new file mode 100644
index 0000000..430a933
Binary files /dev/null and b/13_eigenvalues/output/00011.png differ
diff --git a/13_eigenvalues/output/00012.png b/13_eigenvalues/output/00012.png
new file mode 100644
index 0000000..c708dd3
Binary files /dev/null and b/13_eigenvalues/output/00012.png differ
diff --git a/13_eigenvalues/output/00013.png b/13_eigenvalues/output/00013.png
new file mode 100644
index 0000000..4a25731
Binary files /dev/null and b/13_eigenvalues/output/00013.png differ
diff --git a/13_eigenvalues/output/00014.png b/13_eigenvalues/output/00014.png
new file mode 100644
index 0000000..ea02a38
Binary files /dev/null and b/13_eigenvalues/output/00014.png differ
diff --git a/13_eigenvalues/output/00015.png b/13_eigenvalues/output/00015.png
new file mode 100644
index 0000000..f0f640e
Binary files /dev/null and b/13_eigenvalues/output/00015.png differ
diff --git a/13_eigenvalues/output/00016.png b/13_eigenvalues/output/00016.png
new file mode 100644
index 0000000..44f338f
Binary files /dev/null and b/13_eigenvalues/output/00016.png differ
diff --git a/13_eigenvalues/output/00017.png b/13_eigenvalues/output/00017.png
new file mode 100644
index 0000000..056a85b
Binary files /dev/null and b/13_eigenvalues/output/00017.png differ
diff --git a/13_eigenvalues/output/00018.png b/13_eigenvalues/output/00018.png
new file mode 100644
index 0000000..2bab623
Binary files /dev/null and b/13_eigenvalues/output/00018.png differ
diff --git a/13_eigenvalues/output/00019.png b/13_eigenvalues/output/00019.png
new file mode 100644
index 0000000..9dac5a7
Binary files /dev/null and b/13_eigenvalues/output/00019.png differ
diff --git a/13_eigenvalues/output/00020.png b/13_eigenvalues/output/00020.png
new file mode 100644
index 0000000..1d4ab40
Binary files /dev/null and b/13_eigenvalues/output/00020.png differ
diff --git a/13_eigenvalues/output/00021.png b/13_eigenvalues/output/00021.png
new file mode 100644
index 0000000..66eb44d
Binary files /dev/null and b/13_eigenvalues/output/00021.png differ
diff --git a/13_eigenvalues/output/00022.png b/13_eigenvalues/output/00022.png
new file mode 100644
index 0000000..fa6fdc8
Binary files /dev/null and b/13_eigenvalues/output/00022.png differ
diff --git a/13_eigenvalues/output/00023.png b/13_eigenvalues/output/00023.png
new file mode 100644
index 0000000..aa755e1
Binary files /dev/null and b/13_eigenvalues/output/00023.png differ
diff --git a/13_eigenvalues/output/00024.png b/13_eigenvalues/output/00024.png
new file mode 100644
index 0000000..4e64252
Binary files /dev/null and b/13_eigenvalues/output/00024.png differ
diff --git a/13_eigenvalues/output/00025.png b/13_eigenvalues/output/00025.png
new file mode 100644
index 0000000..a2fc1b3
Binary files /dev/null and b/13_eigenvalues/output/00025.png differ
diff --git a/13_eigenvalues/output/00026.png b/13_eigenvalues/output/00026.png
new file mode 100644
index 0000000..041a839
Binary files /dev/null and b/13_eigenvalues/output/00026.png differ
diff --git a/13_eigenvalues/output/00027.png b/13_eigenvalues/output/00027.png
new file mode 100644
index 0000000..41a54ff
Binary files /dev/null and b/13_eigenvalues/output/00027.png differ
diff --git a/13_eigenvalues/output/00028.png b/13_eigenvalues/output/00028.png
new file mode 100644
index 0000000..b6cc8c8
Binary files /dev/null and b/13_eigenvalues/output/00028.png differ
diff --git a/13_eigenvalues/output/00029.png b/13_eigenvalues/output/00029.png
new file mode 100644
index 0000000..4c665b8
Binary files /dev/null and b/13_eigenvalues/output/00029.png differ
diff --git a/13_eigenvalues/output/00030.png b/13_eigenvalues/output/00030.png
new file mode 100644
index 0000000..71801e9
Binary files /dev/null and b/13_eigenvalues/output/00030.png differ
diff --git a/13_eigenvalues/output/00031.png b/13_eigenvalues/output/00031.png
new file mode 100644
index 0000000..a30c5b9
Binary files /dev/null and b/13_eigenvalues/output/00031.png differ
diff --git a/13_eigenvalues/output/00032.png b/13_eigenvalues/output/00032.png
new file mode 100644
index 0000000..8ef9b30
Binary files /dev/null and b/13_eigenvalues/output/00032.png differ
diff --git a/13_eigenvalues/output/00033.png b/13_eigenvalues/output/00033.png
new file mode 100644
index 0000000..79f8e85
Binary files /dev/null and b/13_eigenvalues/output/00033.png differ
diff --git a/13_eigenvalues/output/00034.png b/13_eigenvalues/output/00034.png
new file mode 100644
index 0000000..075d85e
Binary files /dev/null and b/13_eigenvalues/output/00034.png differ
diff --git a/13_eigenvalues/output/00035.png b/13_eigenvalues/output/00035.png
new file mode 100644
index 0000000..bc73c87
Binary files /dev/null and b/13_eigenvalues/output/00035.png differ
diff --git a/13_eigenvalues/output/00036.png b/13_eigenvalues/output/00036.png
new file mode 100644
index 0000000..48b5fae
Binary files /dev/null and b/13_eigenvalues/output/00036.png differ
diff --git a/13_eigenvalues/output/00037.png b/13_eigenvalues/output/00037.png
new file mode 100644
index 0000000..7107cfc
Binary files /dev/null and b/13_eigenvalues/output/00037.png differ
diff --git a/13_eigenvalues/output/00038.png b/13_eigenvalues/output/00038.png
new file mode 100644
index 0000000..c9c09b2
Binary files /dev/null and b/13_eigenvalues/output/00038.png differ
diff --git a/13_eigenvalues/output/00039.png b/13_eigenvalues/output/00039.png
new file mode 100644
index 0000000..8755cf3
Binary files /dev/null and b/13_eigenvalues/output/00039.png differ
diff --git a/13_eigenvalues/output/00040.png b/13_eigenvalues/output/00040.png
new file mode 100644
index 0000000..9503926
Binary files /dev/null and b/13_eigenvalues/output/00040.png differ
diff --git a/13_eigenvalues/output/00041.png b/13_eigenvalues/output/00041.png
new file mode 100644
index 0000000..4311d25
Binary files /dev/null and b/13_eigenvalues/output/00041.png differ
diff --git a/13_eigenvalues/output/00042.png b/13_eigenvalues/output/00042.png
new file mode 100644
index 0000000..3ce55af
Binary files /dev/null and b/13_eigenvalues/output/00042.png differ
diff --git a/13_eigenvalues/output/00043.png b/13_eigenvalues/output/00043.png
new file mode 100644
index 0000000..8c619cf
Binary files /dev/null and b/13_eigenvalues/output/00043.png differ
diff --git a/13_eigenvalues/output/00044.png b/13_eigenvalues/output/00044.png
new file mode 100644
index 0000000..2d91fd9
Binary files /dev/null and b/13_eigenvalues/output/00044.png differ
diff --git a/13_eigenvalues/output/00045.png b/13_eigenvalues/output/00045.png
new file mode 100644
index 0000000..92e9720
Binary files /dev/null and b/13_eigenvalues/output/00045.png differ
diff --git a/13_eigenvalues/output/00046.png b/13_eigenvalues/output/00046.png
new file mode 100644
index 0000000..62dfa75
Binary files /dev/null and b/13_eigenvalues/output/00046.png differ
diff --git a/13_eigenvalues/output/00047.png b/13_eigenvalues/output/00047.png
new file mode 100644
index 0000000..b62ac38
Binary files /dev/null and b/13_eigenvalues/output/00047.png differ
diff --git a/13_eigenvalues/output/00048.png b/13_eigenvalues/output/00048.png
new file mode 100644
index 0000000..0b90065
Binary files /dev/null and b/13_eigenvalues/output/00048.png differ
diff --git a/13_eigenvalues/output/00049.png b/13_eigenvalues/output/00049.png
new file mode 100644
index 0000000..f7510c0
Binary files /dev/null and b/13_eigenvalues/output/00049.png differ
diff --git a/13_eigenvalues/output/00050.png b/13_eigenvalues/output/00050.png
new file mode 100644
index 0000000..2c74c37
Binary files /dev/null and b/13_eigenvalues/output/00050.png differ
diff --git a/13_eigenvalues/output/00051.png b/13_eigenvalues/output/00051.png
new file mode 100644
index 0000000..daabdee
Binary files /dev/null and b/13_eigenvalues/output/00051.png differ
diff --git a/13_eigenvalues/output/00052.png b/13_eigenvalues/output/00052.png
new file mode 100644
index 0000000..15fb10c
Binary files /dev/null and b/13_eigenvalues/output/00052.png differ
diff --git a/13_eigenvalues/output/00053.png b/13_eigenvalues/output/00053.png
new file mode 100644
index 0000000..07fb58c
Binary files /dev/null and b/13_eigenvalues/output/00053.png differ
diff --git a/13_eigenvalues/output/00054.png b/13_eigenvalues/output/00054.png
new file mode 100644
index 0000000..7bc6b82
Binary files /dev/null and b/13_eigenvalues/output/00054.png differ
diff --git a/13_eigenvalues/output/00055.png b/13_eigenvalues/output/00055.png
new file mode 100644
index 0000000..52df944
Binary files /dev/null and b/13_eigenvalues/output/00055.png differ
diff --git a/13_eigenvalues/output/00056.png b/13_eigenvalues/output/00056.png
new file mode 100644
index 0000000..1677817
Binary files /dev/null and b/13_eigenvalues/output/00056.png differ
diff --git a/13_eigenvalues/output/00057.png b/13_eigenvalues/output/00057.png
new file mode 100644
index 0000000..c4032a1
Binary files /dev/null and b/13_eigenvalues/output/00057.png differ
diff --git a/13_eigenvalues/output/00058.png b/13_eigenvalues/output/00058.png
new file mode 100644
index 0000000..3afae58
Binary files /dev/null and b/13_eigenvalues/output/00058.png differ
diff --git a/13_eigenvalues/output/00059.png b/13_eigenvalues/output/00059.png
new file mode 100644
index 0000000..9e18362
Binary files /dev/null and b/13_eigenvalues/output/00059.png differ
diff --git a/13_eigenvalues/output/00060.png b/13_eigenvalues/output/00060.png
new file mode 100644
index 0000000..47af589
Binary files /dev/null and b/13_eigenvalues/output/00060.png differ
diff --git a/13_eigenvalues/output/00061.png b/13_eigenvalues/output/00061.png
new file mode 100644
index 0000000..b12af4a
Binary files /dev/null and b/13_eigenvalues/output/00061.png differ
diff --git a/13_eigenvalues/output/00062.png b/13_eigenvalues/output/00062.png
new file mode 100644
index 0000000..0999e97
Binary files /dev/null and b/13_eigenvalues/output/00062.png differ
diff --git a/13_eigenvalues/output/00063.png b/13_eigenvalues/output/00063.png
new file mode 100644
index 0000000..5287e29
Binary files /dev/null and b/13_eigenvalues/output/00063.png differ
diff --git a/13_eigenvalues/output/00064.png b/13_eigenvalues/output/00064.png
new file mode 100644
index 0000000..bbad872
Binary files /dev/null and b/13_eigenvalues/output/00064.png differ
diff --git a/13_eigenvalues/output/00065.png b/13_eigenvalues/output/00065.png
new file mode 100644
index 0000000..61d2fe1
Binary files /dev/null and b/13_eigenvalues/output/00065.png differ
diff --git a/13_eigenvalues/output/00066.png b/13_eigenvalues/output/00066.png
new file mode 100644
index 0000000..a1b9eeb
Binary files /dev/null and b/13_eigenvalues/output/00066.png differ
diff --git a/13_eigenvalues/output/00067.png b/13_eigenvalues/output/00067.png
new file mode 100644
index 0000000..67ccce3
Binary files /dev/null and b/13_eigenvalues/output/00067.png differ
diff --git a/13_eigenvalues/output/00068.png b/13_eigenvalues/output/00068.png
new file mode 100644
index 0000000..50c4da5
Binary files /dev/null and b/13_eigenvalues/output/00068.png differ
diff --git a/13_eigenvalues/output/00069.png b/13_eigenvalues/output/00069.png
new file mode 100644
index 0000000..790f009
Binary files /dev/null and b/13_eigenvalues/output/00069.png differ
diff --git a/13_eigenvalues/output/00070.png b/13_eigenvalues/output/00070.png
new file mode 100644
index 0000000..350e64a
Binary files /dev/null and b/13_eigenvalues/output/00070.png differ
diff --git a/13_eigenvalues/output/00071.png b/13_eigenvalues/output/00071.png
new file mode 100644
index 0000000..d4d54b0
Binary files /dev/null and b/13_eigenvalues/output/00071.png differ
diff --git a/13_eigenvalues/output/00072.png b/13_eigenvalues/output/00072.png
new file mode 100644
index 0000000..146d82f
Binary files /dev/null and b/13_eigenvalues/output/00072.png differ
diff --git a/13_eigenvalues/output/00073.png b/13_eigenvalues/output/00073.png
new file mode 100644
index 0000000..cee2012
Binary files /dev/null and b/13_eigenvalues/output/00073.png differ
diff --git a/13_eigenvalues/output/00074.png b/13_eigenvalues/output/00074.png
new file mode 100644
index 0000000..b762ee1
Binary files /dev/null and b/13_eigenvalues/output/00074.png differ
diff --git a/13_eigenvalues/output/00075.png b/13_eigenvalues/output/00075.png
new file mode 100644
index 0000000..59f5ee9
Binary files /dev/null and b/13_eigenvalues/output/00075.png differ
diff --git a/13_eigenvalues/output/00076.png b/13_eigenvalues/output/00076.png
new file mode 100644
index 0000000..59c0066
Binary files /dev/null and b/13_eigenvalues/output/00076.png differ
diff --git a/13_eigenvalues/output/00077.png b/13_eigenvalues/output/00077.png
new file mode 100644
index 0000000..6048a1c
Binary files /dev/null and b/13_eigenvalues/output/00077.png differ
diff --git a/13_eigenvalues/output/00078.png b/13_eigenvalues/output/00078.png
new file mode 100644
index 0000000..4473c93
Binary files /dev/null and b/13_eigenvalues/output/00078.png differ
diff --git a/13_eigenvalues/output/00079.png b/13_eigenvalues/output/00079.png
new file mode 100644
index 0000000..7a6335a
Binary files /dev/null and b/13_eigenvalues/output/00079.png differ
diff --git a/13_eigenvalues/output/00080.png b/13_eigenvalues/output/00080.png
new file mode 100644
index 0000000..821076d
Binary files /dev/null and b/13_eigenvalues/output/00080.png differ
diff --git a/13_eigenvalues/output/00081.png b/13_eigenvalues/output/00081.png
new file mode 100644
index 0000000..d36a35c
Binary files /dev/null and b/13_eigenvalues/output/00081.png differ
diff --git a/13_eigenvalues/output/00082.png b/13_eigenvalues/output/00082.png
new file mode 100644
index 0000000..bde9280
Binary files /dev/null and b/13_eigenvalues/output/00082.png differ
diff --git a/13_eigenvalues/output/00083.png b/13_eigenvalues/output/00083.png
new file mode 100644
index 0000000..1691529
Binary files /dev/null and b/13_eigenvalues/output/00083.png differ
diff --git a/13_eigenvalues/output/00084.png b/13_eigenvalues/output/00084.png
new file mode 100644
index 0000000..5adfef2
Binary files /dev/null and b/13_eigenvalues/output/00084.png differ
diff --git a/13_eigenvalues/output/00085.png b/13_eigenvalues/output/00085.png
new file mode 100644
index 0000000..02067f8
Binary files /dev/null and b/13_eigenvalues/output/00085.png differ
diff --git a/13_eigenvalues/output/00086.png b/13_eigenvalues/output/00086.png
new file mode 100644
index 0000000..2dfc5ac
Binary files /dev/null and b/13_eigenvalues/output/00086.png differ
diff --git a/13_eigenvalues/output/00087.png b/13_eigenvalues/output/00087.png
new file mode 100644
index 0000000..a955fd7
Binary files /dev/null and b/13_eigenvalues/output/00087.png differ
diff --git a/13_eigenvalues/output/00088.png b/13_eigenvalues/output/00088.png
new file mode 100644
index 0000000..7b38164
Binary files /dev/null and b/13_eigenvalues/output/00088.png differ
diff --git a/13_eigenvalues/output/00089.png b/13_eigenvalues/output/00089.png
new file mode 100644
index 0000000..2c75b0f
Binary files /dev/null and b/13_eigenvalues/output/00089.png differ
diff --git a/13_eigenvalues/output/00090.png b/13_eigenvalues/output/00090.png
new file mode 100644
index 0000000..e4489e6
Binary files /dev/null and b/13_eigenvalues/output/00090.png differ
diff --git a/13_eigenvalues/output/00091.png b/13_eigenvalues/output/00091.png
new file mode 100644
index 0000000..af5d877
Binary files /dev/null and b/13_eigenvalues/output/00091.png differ
diff --git a/13_eigenvalues/output/00092.png b/13_eigenvalues/output/00092.png
new file mode 100644
index 0000000..ab960d1
Binary files /dev/null and b/13_eigenvalues/output/00092.png differ
diff --git a/13_eigenvalues/output/00093.png b/13_eigenvalues/output/00093.png
new file mode 100644
index 0000000..81c02c0
Binary files /dev/null and b/13_eigenvalues/output/00093.png differ
diff --git a/13_eigenvalues/output/00094.png b/13_eigenvalues/output/00094.png
new file mode 100644
index 0000000..850739f
Binary files /dev/null and b/13_eigenvalues/output/00094.png differ
diff --git a/13_eigenvalues/output/00095.png b/13_eigenvalues/output/00095.png
new file mode 100644
index 0000000..41d354d
Binary files /dev/null and b/13_eigenvalues/output/00095.png differ
diff --git a/13_eigenvalues/output/00096.png b/13_eigenvalues/output/00096.png
new file mode 100644
index 0000000..7149da5
Binary files /dev/null and b/13_eigenvalues/output/00096.png differ
diff --git a/13_eigenvalues/output/00097.png b/13_eigenvalues/output/00097.png
new file mode 100644
index 0000000..a7bc3f4
Binary files /dev/null and b/13_eigenvalues/output/00097.png differ
diff --git a/13_eigenvalues/output/00098.png b/13_eigenvalues/output/00098.png
new file mode 100644
index 0000000..119aaef
Binary files /dev/null and b/13_eigenvalues/output/00098.png differ
diff --git a/13_eigenvalues/output/00099.png b/13_eigenvalues/output/00099.png
new file mode 100644
index 0000000..f6f7b8a
Binary files /dev/null and b/13_eigenvalues/output/00099.png differ
diff --git a/13_eigenvalues/output/00100.png b/13_eigenvalues/output/00100.png
new file mode 100644
index 0000000..8d37137
Binary files /dev/null and b/13_eigenvalues/output/00100.png differ
diff --git a/13_eigenvalues/springs_masses.png b/13_eigenvalues/springs_masses.png
new file mode 100644
index 0000000..9a01154
Binary files /dev/null and b/13_eigenvalues/springs_masses.png differ
diff --git a/13_eigenvalues/stress.svg b/13_eigenvalues/stress.svg
new file mode 100644
index 0000000..dbbd0b7
--- /dev/null
+++ b/13_eigenvalues/stress.svg
@@ -0,0 +1,297 @@
+
+
+
+
diff --git a/13_eigenvalues/stress_soln.svg b/13_eigenvalues/stress_soln.svg
new file mode 100644
index 0000000..5366d11
--- /dev/null
+++ b/13_eigenvalues/stress_soln.svg
@@ -0,0 +1,297 @@
+
+
+
+
diff --git "a/13_eigenvalues/\360&\341\001" "b/13_eigenvalues/\360&\341\001"
new file mode 100644
index 0000000..e69de29