diff --git a/11_LU-decomposition/mass_springs.svg b/11_LU-decomposition/mass_springs.svg
index c6edff5..abfb4ed 100644
--- a/11_LU-decomposition/mass_springs.svg
+++ b/11_LU-decomposition/mass_springs.svg
@@ -15,7 +15,10 @@
id="svg7576"
version="1.1"
inkscape:version="0.91 r13725"
- sodipodi:docname="mass_springs.svg">
+ sodipodi:docname="mass_springs.svg"
+ inkscape:export-filename="/home/ryan/me3255git/HW4/mass_springs.png"
+ inkscape:export-xdpi="120"
+ inkscape:export-ydpi="120">
diff --git a/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb b/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb
index cc25231..932c3cd 100644
--- a/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb
+++ b/12_matrix_condition/.ipynb_checkpoints/12_matrix_condition-checkpoint.ipynb
@@ -353,4888 +353,6 @@
"\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": {
diff --git a/12_matrix_condition/12_matrix_condition.ipynb b/12_matrix_condition/12_matrix_condition.ipynb
index cc25231..932c3cd 100644
--- a/12_matrix_condition/12_matrix_condition.ipynb
+++ b/12_matrix_condition/12_matrix_condition.ipynb
@@ -353,4888 +353,6 @@
"\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": {
diff --git a/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb b/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb
index f6b2b95..13a92e8 100644
--- a/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb
+++ b/13_eigenvalues/.ipynb_checkpoints/13_eigenvalues-checkpoint.ipynb
@@ -4,7 +4,10 @@
"cell_type": "code",
"execution_count": 1,
"metadata": {
- "collapsed": true
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
},
"outputs": [],
"source": [
@@ -15,7 +18,10 @@
"cell_type": "code",
"execution_count": 2,
"metadata": {
- "collapsed": true
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
},
"outputs": [],
"source": [
@@ -24,22 +30,53 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
- "# Eigenvalues\n",
- "\n",
- "Eigenvalues and eigen vectors are the solution to the set of equations where\n",
+ "# Eigenvalues"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Eigenvalues and eigenvectors 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",
+ "$A-I \\lambda=0$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"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",
+ "These problems are seen in a number of engineering practices:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
"1. Determining vibrational modes in structural devices\n",
"\n",
"2. Material science - vibrational modes of crystal lattices (phonons)\n",
@@ -48,8 +85,17 @@
"\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",
+ "5. Solid mechanics, principle stresses and principle stress directions are eigenvalues and eigenvectors"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
"One way of determining the eigenvalues is taking the determinant:\n",
"\n",
"$|A-\\lambda I|=0$\n",
@@ -59,7 +105,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"Take, A\n",
"\n",
@@ -76,7 +126,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"$(2-\\lambda)(3-\\lambda)(4-\\lambda)-(4-\\lambda)-(2-\\lambda)=0$\n",
"\n",
@@ -98,7 +152,10 @@
"cell_type": "code",
"execution_count": 3,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -122,7 +179,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Applications of Eigenvalue analysis\n",
"\n",
@@ -136,8 +197,17 @@
"\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",
+ "$m_{2}\\frac{d^{2}x_{2}}{dt^{2}}=-k(x_{2}-x_{1})-kx_{2}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"we know that $x_{i}(t)\\propto sin(\\omega t)$ so we can substitute\n",
"\n",
"$x_{i}=X_{i}sin(\\omega t)$\n",
@@ -151,7 +221,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"now, \n",
"\n",
@@ -174,7 +248,10 @@
"cell_type": "code",
"execution_count": 4,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -214,9 +291,12 @@
},
{
"cell_type": "code",
- "execution_count": 5,
+ "execution_count": 4,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -256,7 +336,10 @@
"cell_type": "code",
"execution_count": 6,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -277,24 +360,52 @@
}
],
"source": [
- "K*v(:,2)\n",
+ "K*v(:,2) % K*x=lambda*x\n",
"e(2,2)*v(:,2)"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"### $m_{1}=m_{2}=$40 kg and $k_{1}=k_{2}=k_{3}=$200 N/m\n",
- "![animation](./eig_200_40.gif)\n",
- "\n",
+ "![animation](./eig_200_40.gif)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"### $m_{1}=$40 kg, $m_{2}=$50 kg, $k_{1}=$200 N/m, and $k_{2}=k_{3}=$100 N/m\n",
"![animation](./eig_200_100_40_50.gif)"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Thanks"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Solid Mechanics\n",
"\n",
@@ -567,6 +678,7 @@
}
],
"metadata": {
+ "celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Octave",
"language": "octave",
diff --git a/13_eigenvalues/13_eigenvalues.ipynb b/13_eigenvalues/13_eigenvalues.ipynb
index f6b2b95..13a92e8 100644
--- a/13_eigenvalues/13_eigenvalues.ipynb
+++ b/13_eigenvalues/13_eigenvalues.ipynb
@@ -4,7 +4,10 @@
"cell_type": "code",
"execution_count": 1,
"metadata": {
- "collapsed": true
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
},
"outputs": [],
"source": [
@@ -15,7 +18,10 @@
"cell_type": "code",
"execution_count": 2,
"metadata": {
- "collapsed": true
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
},
"outputs": [],
"source": [
@@ -24,22 +30,53 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
- "# Eigenvalues\n",
- "\n",
- "Eigenvalues and eigen vectors are the solution to the set of equations where\n",
+ "# Eigenvalues"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Eigenvalues and eigenvectors 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",
+ "$A-I \\lambda=0$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"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",
+ "These problems are seen in a number of engineering practices:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
"1. Determining vibrational modes in structural devices\n",
"\n",
"2. Material science - vibrational modes of crystal lattices (phonons)\n",
@@ -48,8 +85,17 @@
"\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",
+ "5. Solid mechanics, principle stresses and principle stress directions are eigenvalues and eigenvectors"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
"One way of determining the eigenvalues is taking the determinant:\n",
"\n",
"$|A-\\lambda I|=0$\n",
@@ -59,7 +105,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"Take, A\n",
"\n",
@@ -76,7 +126,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"$(2-\\lambda)(3-\\lambda)(4-\\lambda)-(4-\\lambda)-(2-\\lambda)=0$\n",
"\n",
@@ -98,7 +152,10 @@
"cell_type": "code",
"execution_count": 3,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -122,7 +179,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Applications of Eigenvalue analysis\n",
"\n",
@@ -136,8 +197,17 @@
"\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",
+ "$m_{2}\\frac{d^{2}x_{2}}{dt^{2}}=-k(x_{2}-x_{1})-kx_{2}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"we know that $x_{i}(t)\\propto sin(\\omega t)$ so we can substitute\n",
"\n",
"$x_{i}=X_{i}sin(\\omega t)$\n",
@@ -151,7 +221,11 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"now, \n",
"\n",
@@ -174,7 +248,10 @@
"cell_type": "code",
"execution_count": 4,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -214,9 +291,12 @@
},
{
"cell_type": "code",
- "execution_count": 5,
+ "execution_count": 4,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -256,7 +336,10 @@
"cell_type": "code",
"execution_count": 6,
"metadata": {
- "collapsed": false
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
},
"outputs": [
{
@@ -277,24 +360,52 @@
}
],
"source": [
- "K*v(:,2)\n",
+ "K*v(:,2) % K*x=lambda*x\n",
"e(2,2)*v(:,2)"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"### $m_{1}=m_{2}=$40 kg and $k_{1}=k_{2}=k_{3}=$200 N/m\n",
- "![animation](./eig_200_40.gif)\n",
- "\n",
+ "![animation](./eig_200_40.gif)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
"### $m_{1}=$40 kg, $m_{2}=$50 kg, $k_{1}=$200 N/m, and $k_{2}=k_{3}=$100 N/m\n",
"![animation](./eig_200_100_40_50.gif)"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Thanks"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Solid Mechanics\n",
"\n",
@@ -567,6 +678,7 @@
}
],
"metadata": {
+ "celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Octave",
"language": "octave",
diff --git a/13_eigenvalues/octave-workspace b/13_eigenvalues/octave-workspace
index 8c437bb..75b054b 100644
Binary files a/13_eigenvalues/octave-workspace and b/13_eigenvalues/octave-workspace differ