diff --git a/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb b/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb index fa0f10b..bcd512f 100644 --- a/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb +++ b/08_optimization/.ipynb_checkpoints/lecture_08-checkpoint.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 10, + "execution_count": 20, "metadata": { "collapsed": false, "slideshow": { @@ -57,8 +57,17 @@ "source": [ "The Lennard-Jones potential is commonly used to model interatomic bonding. \n", "\n", - "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$\n", - "\n", + "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ "Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. \n", "\n", "![TEM image of Gold chain](au_chain.jpg)\n", @@ -88,18 +97,27 @@ } }, "source": [ - "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n", - "\n", + "### First, let's find the minimum energy $\\min(E_{LJ}(x))$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ "## Brute force" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 39, "metadata": { "collapsed": false, "slideshow": { - "slide_type": "subslide" + "slide_type": "fragment" } }, "outputs": [ @@ -240,7 +258,7 @@ "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\n", @@ -263,6 +281,10 @@ } ], "source": [ + "function E_LJ =lennard_jones(x,sigma,epsilon)\n", + " E_LJ = 4*epsilon*((sigma./x).^12-(sigma./x).^6);\n", + "end\n", + "\n", "setdefaults\n", "epsilon = 0.039; % kcal/mol\n", "sigma = 2.934; % Angstrom\n", @@ -277,7 +299,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 40, "metadata": { "collapsed": false, "slideshow": { @@ -322,202 +344,14 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 41, "metadata": { - "collapsed": false, + "collapsed": true, "slideshow": { "slide_type": "subslide" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "current_min = -0.019959\r\n" - ] - }, - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "Gnuplot\n", - "Produced by GNUPLOT 5.0 patchlevel 3 \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t \n", - "\t \n", - "\t\n", - "\t\n", - "\t \n", - "\t \n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\t\n", - "\t\t\n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\t\t\n", - "\t\t-0.04\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t-0.02\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.02\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.04\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.06\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.08\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t2.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t3\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t3.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t4\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t4.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t5.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t6\n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\t\n", - "\n", - "\n", - "\tgnuplot_plot_1a\n", - "\n", - "\n", - "\n", - "\t\n", - "\t\n", - "\tgnuplot_plot_2a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_3a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_4a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_5a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_6a\n", - "\n", - "\t\n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "% define Au atomic potential\n", "epsilon = 0.039; % kcal/mol\n", @@ -535,31 +369,15 @@ "\n", "x1=x_l+d; % define point 1\n", "x2=x_u-d; % define point 2\n", - "\n", - "\n", "% evaluate Au_x(x1) and Au_x(x2)\n", "\n", "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", - "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", - "hold on;\n", - "\n", - "if f2gnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_5a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_6a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -757,17 +574,9 @@ } ], "source": [ - "% Iteration #2\n", - "d=(phi-1)*(x_u-x_l);\n", - "\n", - "x1=x_l+d; % define point 1\n", - "x2=x_u-d; % define point 2\n", - "\n", - "% evaluate Au_x(x1) and Au_x(x2)\n", "\n", - "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", - "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", + "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',...\n", + "x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", "hold on;\n", "\n", "if f2gnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_5a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_6a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -983,19 +810,8 @@ } ], "source": [ - "% Iteration #3\n", - "d=(phi-1)*(x_u-x_l);\n", - "\n", - "x1=x_l+d; % define point 1\n", - "x2=x_u-d; % define point 2\n", - "\n", - "% evaluate Au_x(x1) and Au_x(x2)\n", - "\n", - "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", "hold on;\n", - "\n", "if f2gnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_5a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_6a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -1209,19 +1048,8 @@ } ], "source": [ - "% Iteration #3\n", - "d=(phi-1)*(x_u-x_l);\n", - "\n", - "x1=x_l+d; % define point 1\n", - "x2=x_u-d; % define point 2\n", - "\n", - "% evaluate Au_x(x1) and Au_x(x2)\n", - "\n", - "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", "hold on;\n", - "\n", "if f2= maxit,break,end\n", + "end\n", + "x=xopt;fx=f(xopt,varargin{:});" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Minimization with goldmin" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "xopt = 3.2933\n", + "fmin = -0.039000\n", + "ea = 6.3350e-04\n", + "iters = 23\n" + ] + } + ], + "source": [ + "[xopt,fmin,ea,iters]=goldmin(Au_x,2.5,6,0.001)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, { "cell_type": "markdown", "metadata": { @@ -1924,7 +1867,9 @@ } }, "source": [ - "Parabolic interpolation does not converge in many scenarios even though it it a bracketing method. Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)" + "Parabolic interpolation does not always converge despite being a bracketing method. \n", + "\n", + "Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)" ] }, { @@ -1946,7 +1891,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 5, "metadata": { "collapsed": false, "slideshow": { @@ -1958,11 +1903,6 @@ "name": "stdout", "output_type": "stream", "text": [ - "warning: Matlab-style short-circuit operation performed for operator |\n", - "warning: called from\n", - " goldmin at line 17 column 1\n", - "warning: Matlab-style short-circuit operation performed for operator |\n", - "warning: Matlab-style short-circuit operation performed for operator |\n", "xmin = 0.32933\n", "Emin = -2.7096e-04\n" ] @@ -2099,7 +2039,7 @@ "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\n", @@ -2130,11 +2070,8 @@ "sigma = 2.934; % Angstrom\n", "sigma = sigma*0.10; % nm/Angstrom\n", "x=linspace(2.8,6,200)*0.10; % bond length in um\n", - "\n", "Ex = lennard_jones(x,sigma,epsilon);\n", "\n", - "%[Emin,imin]=min(Ex);\n", - "\n", "[xmin,Emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n", "\n", "plot(x,Ex,xmin,Emin,'o')\n", diff --git a/08_optimization/lecture_08.ipynb b/08_optimization/lecture_08.ipynb index fa0f10b..bcd512f 100644 --- a/08_optimization/lecture_08.ipynb +++ b/08_optimization/lecture_08.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 10, + "execution_count": 20, "metadata": { "collapsed": false, "slideshow": { @@ -57,8 +57,17 @@ "source": [ "The Lennard-Jones potential is commonly used to model interatomic bonding. \n", "\n", - "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$\n", - "\n", + "$E_{LJ}(x)=4\\epsilon \\left(\\left(\\frac{\\sigma}{x}\\right)^{12}-\\left(\\frac{\\sigma}{x}\\right)^{6}\\right)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ "Considering a 1-D gold chain, we can calculate the bond length, $x_{b}$, with no force applied to the chain and even for tension, F. This will allow us to calculate the nonlinear spring constant of a 1-D gold chain. \n", "\n", "![TEM image of Gold chain](au_chain.jpg)\n", @@ -88,18 +97,27 @@ } }, "source": [ - "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n", - "\n", + "### First, let's find the minimum energy $\\min(E_{LJ}(x))$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ "## Brute force" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 39, "metadata": { "collapsed": false, "slideshow": { - "slide_type": "subslide" + "slide_type": "fragment" } }, "outputs": [ @@ -240,7 +258,7 @@ "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\n", @@ -263,6 +281,10 @@ } ], "source": [ + "function E_LJ =lennard_jones(x,sigma,epsilon)\n", + " E_LJ = 4*epsilon*((sigma./x).^12-(sigma./x).^6);\n", + "end\n", + "\n", "setdefaults\n", "epsilon = 0.039; % kcal/mol\n", "sigma = 2.934; % Angstrom\n", @@ -277,7 +299,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 40, "metadata": { "collapsed": false, "slideshow": { @@ -322,202 +344,14 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 41, "metadata": { - "collapsed": false, + "collapsed": true, "slideshow": { "slide_type": "subslide" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "current_min = -0.019959\r\n" - ] - }, - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "Gnuplot\n", - "Produced by GNUPLOT 5.0 patchlevel 3 \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t \n", - "\t \n", - "\t\n", - "\t\n", - "\t \n", - "\t \n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\t\n", - "\t\t\n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\t\t\n", - "\t\t-0.04\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t-0.02\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.02\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.04\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.06\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t0.08\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t2.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t3\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t3.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t4\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t4.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t5.5\n", - "\t\n", - "\n", - "\n", - "\t\t\n", - "\t\t6\n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\t\n", - "\n", - "\n", - "\tgnuplot_plot_1a\n", - "\n", - "\n", - "\n", - "\t\n", - "\t\n", - "\tgnuplot_plot_2a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_3a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_4a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_5a\n", - "\n", - "\t \n", - "\t\n", - "\n", - "\t\n", - "\tgnuplot_plot_6a\n", - "\n", - "\t\n", - "\t\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "% define Au atomic potential\n", "epsilon = 0.039; % kcal/mol\n", @@ -535,31 +369,15 @@ "\n", "x1=x_l+d; % define point 1\n", "x2=x_u-d; % define point 2\n", - "\n", - "\n", "% evaluate Au_x(x1) and Au_x(x2)\n", "\n", "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", - "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", - "hold on;\n", - "\n", - "if f2gnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_5a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_6a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -757,17 +574,9 @@ } ], "source": [ - "% Iteration #2\n", - "d=(phi-1)*(x_u-x_l);\n", - "\n", - "x1=x_l+d; % define point 1\n", - "x2=x_u-d; % define point 2\n", - "\n", - "% evaluate Au_x(x1) and Au_x(x2)\n", "\n", - "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", - "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", + "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',...\n", + "x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", "hold on;\n", "\n", "if f2gnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_5a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_6a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -983,19 +810,8 @@ } ], "source": [ - "% Iteration #3\n", - "d=(phi-1)*(x_u-x_l);\n", - "\n", - "x1=x_l+d; % define point 1\n", - "x2=x_u-d; % define point 2\n", - "\n", - "% evaluate Au_x(x1) and Au_x(x2)\n", - "\n", - "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", "hold on;\n", - "\n", "if f2gnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_5a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_6a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -1209,19 +1048,8 @@ } ], "source": [ - "% Iteration #3\n", - "d=(phi-1)*(x_u-x_l);\n", - "\n", - "x1=x_l+d; % define point 1\n", - "x2=x_u-d; % define point 2\n", - "\n", - "% evaluate Au_x(x1) and Au_x(x2)\n", - "\n", - "f1=Au_x(x1);\n", - "f2=Au_x(x2);\n", "plot(x,Au_x(x),x_l,Au_x(x_l),'ro',x2,f2,'rs',x1,f1,'gs',x_u,Au_x(x_u),'go')\n", "hold on;\n", - "\n", "if f2= maxit,break,end\n", + "end\n", + "x=xopt;fx=f(xopt,varargin{:});" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Minimization with goldmin" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "xopt = 3.2933\n", + "fmin = -0.039000\n", + "ea = 6.3350e-04\n", + "iters = 23\n" + ] + } + ], + "source": [ + "[xopt,fmin,ea,iters]=goldmin(Au_x,2.5,6,0.001)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + }, { "cell_type": "markdown", "metadata": { @@ -1924,7 +1867,9 @@ } }, "source": [ - "Parabolic interpolation does not converge in many scenarios even though it it a bracketing method. Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)" + "Parabolic interpolation does not always converge despite being a bracketing method. \n", + "\n", + "Instead, functions like `fminbnd` in Matlab and Octave use a combination of the two (Golden Ratio and Parabolic)" ] }, { @@ -1946,7 +1891,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 5, "metadata": { "collapsed": false, "slideshow": { @@ -1958,11 +1903,6 @@ "name": "stdout", "output_type": "stream", "text": [ - "warning: Matlab-style short-circuit operation performed for operator |\n", - "warning: called from\n", - " goldmin at line 17 column 1\n", - "warning: Matlab-style short-circuit operation performed for operator |\n", - "warning: Matlab-style short-circuit operation performed for operator |\n", "xmin = 0.32933\n", "Emin = -2.7096e-04\n" ] @@ -2099,7 +2039,7 @@ "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\n", @@ -2130,11 +2070,8 @@ "sigma = 2.934; % Angstrom\n", "sigma = sigma*0.10; % nm/Angstrom\n", "x=linspace(2.8,6,200)*0.10; % bond length in um\n", - "\n", "Ex = lennard_jones(x,sigma,epsilon);\n", "\n", - "%[Emin,imin]=min(Ex);\n", - "\n", "[xmin,Emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n", "\n", "plot(x,Ex,xmin,Emin,'o')\n", diff --git a/08_optimization/octave-workspace b/08_optimization/octave-workspace index 8c437bb..4df051d 100644 Binary files a/08_optimization/octave-workspace and b/08_optimization/octave-workspace differ diff --git a/HW2/README.html b/HW2/README.html new file mode 100644 index 0000000..078c879 --- /dev/null +++ b/HW2/README.html @@ -0,0 +1,77 @@ + + + + + + + + + + + +

Homework #2

+

due 10/6/17 by 11:59pm

+

1. Create a new github repository called ‘02_roots_and_optimization’.

+
    +
  1. Add rcc02007 and zhs15101 as collaborators.

  2. +
  3. submit the clone repository URL to: https://goo.gl/forms/svFKpfiCfLO9Zvfz1

  4. +
+

2. You’re installing a powerline in a residential neighborhood. The lowest point on the cable is 30 m above the ground, but 30 m away is a tree that is 35 m tall. Another engineer informs you that this is a catenary cable problem with the following solution

+
+eq. 1 +

eq. 1

+
+

\(y(x)=\frac{T}{w}\cosh\left(\frac{w}{T}x\right)+y_{0}-\frac{T}{w}\).

+

where y(x) is the height of the cable at a distance, x, from the lowest point, \(y_{0}\), T is the tension in the cable, and w is the weight per unit length of the cable. Your supervisor wants to know which numerical solver to use when they have to install these powerlines in similar places.

+
    +
  1. Use the three solvers falsepos.m, bisect.m, and mod_secant.m to solve for the tension neededi, T, to reach y(30 m)=35 m, with w=10 N/m, and \(y_{0}\)=30 m.

  2. +
  3. Compare the number of iterations that each function needed to reach an accuracy of 0.00001%. Include a table in your README.md with:

    +
    | solver | initial guess(es) | ea | number of iterations|
    +| --- | --- | --- | --- |
    +|falsepos   |  |  |  |
    +|mod_secant |  |  |  |
    +|bisect     |  |  |  |
  4. +
  5. Add a figure to your README that plots the final shape of the powerline (eq2) from x=-10 to 50 m.

  6. +
+

3. The Newton-Raphson method and the modified secant method do not always converge to a solution. One simple example is the function f(x) = (x-1)*exp(-(x-1)^2). The root is at 1, but using the numerical solvers, newtraph.m and mod_secant.m, there are certain initial guesses that do not converge.

+
    +
  1. Calculate the first 5 iterations for the Newton-Raphson method with an initial guess of x_i=3 for f(x)=(x-1)*exp(-(x-1)^2).

  2. +
  3. Add the results to a table in the README.md with:

    +
    ### divergence of Newton-Raphson method
    +
    +| iteration | x_i | approx error |
    +| --- | --- | --- |
    +| 0 | 3 | n/a |
    +| 1 |   |     |
    +| 2 |   |     |
    +| 3 |   |     |
    +| 4 |   |     |
    +| 5 |   |     |
  4. +
  5. Repeat steps a-b for an initial guess of 1.2. (But change the heading from ‘divergence’ to ‘convergence’)

  6. +
+
+Model of Gold chain, from molecular dynamics simulation +

Model of Gold chain, from molecular dynamics simulation

+
+

4. Determine the nonlinear spring constants of a single-atom gold chain. You can assume the gold atoms are aligned in a one dimensional network and the potential energy is described by the Lennard-Jones potential as such

+
+eq3 +

eq3

+
+

\(E_{LJ}(x)=4\epsilon \left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)\).

+

Where x is the distance between atoms in nm, \(\epsilon\)=2.71E-4 aJ, and \(\sigma\)=0.2934 nm. The energy term that must be minimized is

+
+eq4 +

eq4

+
+

\(E_{total}(\Delta x)=E_{LJ}(x_{0}+\Delta x)-F\Delta x\).

+

Where \(x_{0}\) is the distance between atoms with no force applied and dx is the amount each gold atom has moved under a given force, F.

+
    +
  1. Determine x0 when F=0 nN using the golden ratio and parabolic methods. Show your script and output in your README and include your functions

  2. +
  3. Solve for dx is the amount each gold atom has mov for F=0 to 0.0022 nN with 30 steps. *Use the golden ratio solver or the matlab/octave fminsearch

  4. +
  5. create a sum of squares error function sse_of_parabola.m that calculates the sum of squares error between a function \(F(x)=K_{1}\Delta x+1/2K_{2}\Delta x^{2}\) and the Forces used in part B for each dx.

  6. +
  7. Use the fminsearch matlab/octave function to determine k1k2.

  8. +
  9. Plot the force vs calculated dx and the best-fit parabola using k1k2 in part d.

  10. +
+ + diff --git a/HW2/README.md b/HW2/README.md index a3b064d..0b5d703 100644 --- a/HW2/README.md +++ b/HW2/README.md @@ -14,6 +14,8 @@ cable is 30 m above the ground, but 30 m away is a tree that is 35 m tall. Another engineer informs you that this is a catenary cable problem with the following solution + ![eq. 1](./equations/eq1.png) + $y(x)=\frac{T}{w}\cosh\left(\frac{w}{T}x\right)+y_{0}-\frac{T}{w}$. where y(x) is the height of the cable at a distance, x, from the lowest point, $y_{0}$, @@ -37,7 +39,7 @@ engineer informs you that this is a catenary cable problem with the following so c. Add a figure to your README that plots the final shape of the powerline - ($y(x)~vs.~x$) from x=-10 to 50 m. + (![eq2](./equations/eq2.png)) from x=-10 to 50 m. **3\.** The Newton-Raphson method and the modified secant method do not always converge to a solution. One simple example is the function f(x) = (x-1)*exp(-(x-1)^2). The root is at 1, but @@ -71,29 +73,35 @@ guesses that do not converge. the gold atoms are aligned in a one dimensional network and the potential energy is described by the Lennard-Jones potential as such + ![eq3](./equations/eq3.png) + $E_{LJ}(x)=4\epsilon \left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)$. Where x is the distance between atoms in nm, $\epsilon$=2.71E-4 aJ, and $\sigma$=0.2934 nm. The energy term that must be minimized is + ![eq4](./equations/eq4.png) + $E_{total}(\Delta x)=E_{LJ}(x_{0}+\Delta x)-F\Delta x$. - Where $x_{0}$ is the distance between atoms with no force applied and $\Delta x$ is the + Where ![x0](./equations/x0.png) is the distance between atoms with no force applied and + ![dx](./equations/deltax.png) is the amount each gold atom has moved under a given force, F. - a. Determine $x_{0}$ when F=0 nN using the golden ratio and parabolic methods. *Show + a. Determine ![x0](./equations/x0.png) when F=0 nN using the golden ratio and parabolic methods. *Show your script and output in your README and include your functions* - b. Solve for $\Delta x$ for F=0 to 0.0022 nN with 30 steps. *Use the golden ratio + b. Solve for ![dx](./equations/deltax.png) is the + amount each gold atom has mov for F=0 to 0.0022 nN with 30 steps. *Use the golden ratio solver or the matlab/octave `fminsearch` c. create a sum of squares error function `sse_of_parabola.m` that calculates the sum of squares error between a function $F(x)=K_{1}\Delta x+1/2K_{2}\Delta x^{2}$ and the - Forces used in part B for each $\Delta x$. + Forces used in part B for each ![dx](./equations/deltax.png). - d. Use the `fminsearch` matlab/octave function to determine $K_{1}$ and $K_{2}$. + d. Use the `fminsearch` matlab/octave function to determine + ![k1k2](./equations/k1k2.png). - e. Plot the force vs calculated $\Delta x$ and the best-fit parabola using $K_{1}$ and - $K_{2}$ in part d. + e. Plot the force vs calculated ![dx](./equations/deltax.png) and the best-fit parabola using ![k1k2](./equations/k1k2.png) in part d. diff --git a/HW2/equations/README.md b/HW2/equations/README.md new file mode 100644 index 0000000..e69de29 diff --git a/HW2/equations/deltax.png b/HW2/equations/deltax.png new file mode 100644 index 0000000..a27cd7a Binary files /dev/null and b/HW2/equations/deltax.png differ diff --git a/HW2/equations/deltax.tex b/HW2/equations/deltax.tex new file mode 100644 index 0000000..dfba7b4 --- /dev/null +++ b/HW2/equations/deltax.tex @@ -0,0 +1 @@ +$\Delta x$ diff --git a/HW2/equations/eq1.png b/HW2/equations/eq1.png new file mode 100644 index 0000000..fe78ec6 Binary files /dev/null and b/HW2/equations/eq1.png differ diff --git a/HW2/equations/eq1.tex b/HW2/equations/eq1.tex new file mode 100644 index 0000000..3dadcac --- /dev/null +++ b/HW2/equations/eq1.tex @@ -0,0 +1,2 @@ +$y(x)=T/w \cosh\left(\frac{w}{T}x\right)+y_{0}-\frac{T}{w}$. + diff --git a/HW2/equations/eq2.png b/HW2/equations/eq2.png new file mode 100644 index 0000000..2bd0218 Binary files /dev/null and b/HW2/equations/eq2.png differ diff --git a/HW2/equations/eq2.tex b/HW2/equations/eq2.tex new file mode 100644 index 0000000..d6debd4 --- /dev/null +++ b/HW2/equations/eq2.tex @@ -0,0 +1,2 @@ + ($y(x)~vs.~x$) + diff --git a/HW2/equations/eq3.png b/HW2/equations/eq3.png new file mode 100644 index 0000000..121bdf0 Binary files /dev/null and b/HW2/equations/eq3.png differ diff --git a/HW2/equations/eq3.tex b/HW2/equations/eq3.tex new file mode 100644 index 0000000..b807515 --- /dev/null +++ b/HW2/equations/eq3.tex @@ -0,0 +1,3 @@ +$E_{LJ}(x)=4\epsilon +\left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)$. + diff --git a/HW2/equations/eq4.png b/HW2/equations/eq4.png new file mode 100644 index 0000000..2077360 Binary files /dev/null and b/HW2/equations/eq4.png differ diff --git a/HW2/equations/eq4.tex b/HW2/equations/eq4.tex new file mode 100644 index 0000000..26fe469 --- /dev/null +++ b/HW2/equations/eq4.tex @@ -0,0 +1,2 @@ +$E_{total}(\Delta x)=E_{LJ}(x_{0}+\Delta x)-F\Delta x$. + diff --git a/HW2/equations/fx.png b/HW2/equations/fx.png new file mode 100644 index 0000000..bd17982 Binary files /dev/null and b/HW2/equations/fx.png differ diff --git a/HW2/equations/fx.tex b/HW2/equations/fx.tex new file mode 100644 index 0000000..07d410c --- /dev/null +++ b/HW2/equations/fx.tex @@ -0,0 +1,2 @@ +$F(x)=K_{1}\Delta x+1/2K_{2}\Delta x^{2}$ + diff --git a/HW2/equations/k1k2.png b/HW2/equations/k1k2.png new file mode 100644 index 0000000..d7048f7 Binary files /dev/null and b/HW2/equations/k1k2.png differ diff --git a/HW2/equations/k1k2.tex b/HW2/equations/k1k2.tex new file mode 100644 index 0000000..c0d4c95 --- /dev/null +++ b/HW2/equations/k1k2.tex @@ -0,0 +1 @@ +$K_{1}$ and $K_{2}$ diff --git a/HW2/equations/x0.png b/HW2/equations/x0.png new file mode 100644 index 0000000..468b8ff Binary files /dev/null and b/HW2/equations/x0.png differ diff --git a/HW2/equations/x0.tex b/HW2/equations/x0.tex new file mode 100644 index 0000000..037c2d3 --- /dev/null +++ b/HW2/equations/x0.tex @@ -0,0 +1 @@ +$x_{0}$