diff --git a/HW3/README.md b/HW3/README.md index 494e0d5..45214a4 100644 --- a/HW3/README.md +++ b/HW3/README.md @@ -11,8 +11,8 @@ c. Copy your `projectile.m` function into the 'roots_and_optimization' folder. *Disable the plotting routine for the solvers* - d. Use the four solvers `falsepos.m`, `incsearch.m`, `newtraph.m` and `mod_secant.m` - to solve for the angle needed to reach h=1.72 m, with an initial speed of 1.5 m/s. + d. Use the four solvers `falsepos.m`, `bisect.m`, `newtraph.m` and `mod_secant.m` + to solve for the angle needed to reach h=1.72 m, with an initial speed of 15 m/s. e. The `newtraph.m` function needs a derivative, calculate the derivative of your function with respect to theta, `dprojectile_dtheta.m`. This function should @@ -29,7 +29,7 @@ | solver | initial guess(es) | ea | number of iterations| | --- | --- | --- | --- | |falsepos | | | | - |incsearch | | | | + |bisect | | | | |newtraph | | | | |mod_secant | | | | ``` @@ -49,7 +49,7 @@ using the numerical solvers, `newtraph.m` and `mod_secant.m`, there are certain guesses that do not converge. a. Calculate the first 5 iterations for the Newton-Raphson method with an initial - guess of x_i=2. + guess of x_i=2 for f(x)=x*exp(-x^2). b. Add the results to a table in the `README.md` with: @@ -70,5 +70,4 @@ guesses that do not converge. 'divergence' to 'convergence') 3. Commit your changes to your repository. Sync your local repository with github. Then -copy and paste the "clone URL" into the following Google Form [Homework -#3](https://goo.gl/forms/UJBGwp0fQcSxImkq2) +copy and paste the "clone URL" into the following Google Form [Homework 3](https://goo.gl/forms/UJBGwp0fQcSxImkq2) diff --git a/lecture_07/.lecture_07.md.swp b/lecture_07/.lecture_07.md.swp new file mode 100644 index 0000000..0998c86 Binary files /dev/null and b/lecture_07/.lecture_07.md.swp differ diff --git a/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb b/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb new file mode 100644 index 0000000..2fd6442 --- /dev/null +++ b/lecture_08/.ipynb_checkpoints/lecture_08-checkpoint.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_08/Auchain_model.gif b/lecture_08/Auchain_model.gif new file mode 100644 index 0000000..2f8c78c Binary files /dev/null and b/lecture_08/Auchain_model.gif differ diff --git a/lecture_08/au_chain.jpg b/lecture_08/au_chain.jpg new file mode 100644 index 0000000..492fee8 Binary files /dev/null and b/lecture_08/au_chain.jpg differ diff --git a/lecture_08/goldenratio.png b/lecture_08/goldenratio.png new file mode 100644 index 0000000..9493c8c Binary files /dev/null and b/lecture_08/goldenratio.png differ diff --git a/lecture_08/goldmin.m b/lecture_08/goldmin.m new file mode 100644 index 0000000..e5b0ecf --- /dev/null +++ b/lecture_08/goldmin.m @@ -0,0 +1,36 @@ +function [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,varargin) +% goldmin: minimization golden section search +% [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,p1,p2,...): +% uses golden section search to find the minimum of f +% input: +% f = name of function +% xl, xu = lower and upper guesses +% es = desired relative error (default = 0.0001%) +% maxit = maximum allowable iterations (default = 50) +% p1,p2,... = additional parameters used by f +% output: +% x = location of minimum +% fx = minimum function value +% ea = approximate relative error (%) +% iter = number of iterations +if nargin<3,error('at least 3 input arguments required'),end +if nargin<4|isempty(es), es=0.0001;end +if nargin<5|isempty(maxit), maxit=50;end +phi=(1+sqrt(5))/2; +iter=0; +while(1) + d = (phi-1)*(xu - xl); + x1 = xl + d; + x2 = xu - d; + if f(x1,varargin{:}) < f(x2,varargin{:}) + xopt = x1; + xl = x2; + else + xopt = x2; + xu = x1; + end + iter=iter+1; + if xopt~=0, ea = (2 - phi) * abs((xu - xl) / xopt) * 100;end + if ea <= es | iter >= maxit,break,end +end +x=xopt;fx=f(xopt,varargin{:}); \ No newline at end of file diff --git a/lecture_08/lecture_08.ipynb b/lecture_08/lecture_08.ipynb new file mode 100644 index 0000000..d99107a --- /dev/null +++ b/lecture_08/lecture_08.ipynb @@ -0,0 +1,2847 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Optimization\n", + "\n", + "Many problems involve finding a minimum or maximum based on given constraints. Engineers and scientists typically use energy balance equations to find the conditions of minimum energy, but this value may never go to 0 and the actual value of energy may not be of interest. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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", + "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", + "\n", + "Computational Tools to Study and Predict the Long-Term Stability of Nanowires.\n", + "By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva \n", + "\n", + "![Model of Gold chain, from molecular dynamics simulation](Auchain_model.gif)\n", + "\n", + "[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/)\n", + "\n", + "### First, let's find the minimum energy $\\min(E_{LJ}(x))$\n", + "\n", + "## Brute force" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "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", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "setdefaults\n", + "epsilon = 0.039; % kcal/mol\n", + "sigma = 2.934; % Angstrom\n", + "x=linspace(2.8,6,200); % bond length in Angstrom\n", + "\n", + "Ex = lennard_jones(x,sigma,epsilon);\n", + "\n", + "[Emin,imin]=min(Ex);\n", + "\n", + "plot(x,Ex,x(imin),Emin,'o')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 3.2824\n", + "ans = 3.2985\n", + "ans = 3.3146\n" + ] + } + ], + "source": [ + "x(imin-1)\n", + "x(imin)\n", + "x(imin+1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Golden Search Algorithm\n", + "\n", + "We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds.\n", + "\n", + "Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. \n", + "\n", + "![golden ratio](goldenratio.png)\n", + "\n", + "Start with bounds of 2.5 and 6 Angstrom. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "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" + } + ], + "source": [ + "% define Au atomic potential\n", + "epsilon = 0.039; % kcal/mol\n", + "sigma = 2.934; % Angstrom\n", + "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n", + "\n", + "% calculate golden ratio\n", + "phi = 1/2+sqrt(5)/2;\n", + "% set initial limits\n", + "x_l=2.8; \n", + "x_u=6; \n", + "\n", + "% Iteration #1\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", + "\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\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" + } + ], + "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", + "hold on;\n", + "\n", + "if f2\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" + } + ], + "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\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" + } + ], + "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\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.06\n", + "\t\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", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\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", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "% define Au atomic potential\n", + "epsilon = 0.039; % kcal/mol\n", + "sigma = 2.934; % Angstrom\n", + "Au_x= @(x) lennard_jones(x,sigma,epsilon);\n", + "\n", + "% set initial limits\n", + "x_l=2.8; \n", + "x_u=3.5; \n", + "\n", + "% Iteration #1\n", + "x1=x_l;\n", + "x2=mean([x_l,x_u]);\n", + "x3=x_u;\n", + "\n", + "% evaluate Au_x(x1), Au_x(x2) and Au_x(x3)\n", + " \n", + "f1=Au_x(x1);\n", + "f2=Au_x(x2);\n", + "f3=Au_x(x3);\n", + "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n", + "x_fit = linspace(x1,x3,20);\n", + "y_fit = polyval(p,x_fit);\n", + "\n", + "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n", + "hold on\n", + "if f2x2\n", + " plot(x4,f4,'*',[x1,x2],[f1,f2])\n", + " x1=x2;\n", + " f1=f2;\n", + " else\n", + " plot(x4,f4,'*',[x3,x2],[f3,f2])\n", + " x3=x2;\n", + " f3=f2;\n", + " end\n", + " x2=x4; f2=f4;\n", + "else\n", + " error('no minimum in bracket')\n", + "end\n", + "hold off" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "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", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\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", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n", + "x_fit = linspace(x1,x3,20);\n", + "y_fit = polyval(p,x_fit);\n", + "\n", + "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n", + "hold on\n", + "if f2x2\n", + " plot(x4,f4,'*',[x1,x2],[f1,f2])\n", + " x1=x2;\n", + " f1=f2;\n", + " else\n", + " plot(x4,f4,'*',[x3,x2],[f3,f2])\n", + " x3=x2;\n", + " f3=f2;\n", + " end\n", + " x2=x4; f2=f4;\n", + "else\n", + " error('no minimum in bracket')\n", + "end\n", + "hold off\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "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", + "\tgnuplot_plot_3a\n", + "\n", + "\t \n", + "\t\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", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "p = polyfit([x1,x2,x3],[f1,f2,f3],2);\n", + "x_fit = linspace(x1,x3,20);\n", + "y_fit = polyval(p,x_fit);\n", + "\n", + "plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')\n", + "hold on\n", + "if f2x2\n", + " plot(x4,f4,'*',[x1,x2],[f1,f2])\n", + " x1=x2;\n", + " f1=f2;\n", + " else\n", + " plot(x4,f4,'*',[x3,x2],[f3,f2])\n", + " x3=x2;\n", + " f3=f2;\n", + " end\n", + " x2=x4; f2=f4;\n", + "else\n", + " error('no minimum in bracket')\n", + "end\n", + "hold off\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Using the solutions to minimization for the nonlinear spring constant\n", + "\n", + "Now, we have two routines to find minimums of a univariate function (Golden Ratio and Parabolic). Let's use these to solve for the minimum energy given a range of applied forces to the single atom gold chain\n", + "\n", + "$E_{total}(\\Delta x) = E_{LJ}(x_{min}+\\Delta x) - F \\cdot \\Delta x$" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "xmin = 0.32933\n", + "Emin = -2.7096e-04\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.0004\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.0002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0004\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0006\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.7\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tLennard Jones Potential (aJ/atom)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tbond length (nm)\n", + "\t\n", + "\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", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "epsilon = 0.039; % kcal/mol\n", + "epsilon = epsilon*6.9477e-21; % J/atom\n", + "epsilon = epsilon*1e18; % aJ/J\n", + "% final units for epsilon are aJ\n", + "\n", + "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] = fminbnd(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)\n", + "\n", + "plot(x,Ex,xmin,Emin,'o')\n", + "ylabel('Lennard Jones Potential (aJ/atom)')\n", + "xlabel('bond length (nm)')\n", + "\n", + "Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;\n" + ] + }, + { + "cell_type": "code", + "execution_count": 114, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "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", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\n", + "warning: Matlab-style short-circuit operation performed for operator |\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.03\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tForce (nN)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tdx (nm)\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N=50;\n", + "dx = zeros(1,N); % [in nm]\n", + "F_applied=linspace(0,0.0022,N); % [in nN]\n", + "for i=1:N\n", + " optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035);\n", + " dx(i)=optmin;\n", + "end\n", + "\n", + "plot(dx,F_applied)\n", + "xlabel('dx (nm)')\n", + "ylabel('Force (nN)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this function, it is possible to take a derivative and compare the analytical result:" + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "F =\n", + "\n", + "@(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13)\n", + "\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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.03\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.05\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.06\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", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "dx_full=linspace(0,0.06,50);\n", + "F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13)\n", + "plot(dx_full,F(dx_full),dx,F_applied)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Curve-fitting\n", + "Another example is minimizing error in your approximation of a function. If you have data (now we have Force-displacement data) we can fit this to a function, such as:\n", + "\n", + "$F(x) = K_{1}\\Delta x + \\frac{1}{2} K_{2}(\\Delta x)^{2}$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 106, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "function SSE = sse_of_parabola(K,xdata,ydata)\n", + " % calculate the sum of squares error for a parabola given a function, func, and xdata and ydata\n", + " % output is SSE=sum of squares error\n", + " K1=K(1);\n", + " K2=K(2);\n", + " y_function = K1*xdata+1/2*K2*xdata.^2;\n", + " SSE = sum((ydata-y_function).^2);\n", + "end\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 113, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.98 nN/nm^2\n", + "The mininum sum of squares error = 7.35e-08\n" + ] + } + ], + "source": [ + "[K,SSE_min]=fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]);\n", + "fprintf('\\nThe nonlinear spring constants are K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\\n',K)\n", + "fprintf('The mininum sum of squares error = %1.2e',SSE_min)" + ] + }, + { + "cell_type": "code", + "execution_count": 108, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "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\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.001\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.002\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.0025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.005\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.01\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.015\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.02\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.025\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.03\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.035\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.04\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\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", + "\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", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)" + ] + } + ], + "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/lecture_08/lecture_08.md b/lecture_08/lecture_08.md new file mode 100644 index 0000000..8a81a2a --- /dev/null +++ b/lecture_08/lecture_08.md @@ -0,0 +1,626 @@ + + +```octave +%plot --format svg +``` + +# Optimization + +Many problems involve finding a minimum or maximum based on given constraints. Engineers and scientists typically use energy balance equations to find the conditions of minimum energy, but this value may never go to 0 and the actual value of energy may not be of interest. + +The Lennard-Jones potential is commonly used to model interatomic bonding. + +$E_{LJ}(x)=4\epsilon \left(\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right)$ + +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. + +![TEM image of Gold chain](au_chain.jpg) + +Computational Tools to Study and Predict the Long-Term Stability of Nanowires. +By Martin E. Zoloff Michoff, Patricio Vélez, Sergio A. Dassie and Ezequiel P. M. Leiva + +![Model of Gold chain, from molecular dynamics simulation](Auchain_model.gif) + +[Single atom gold chain mechanics](http://www.uam.es/personal_pdi/ciencias/agrait/) + +### First, let's find the minimum energy $\min(E_{LJ}(x))$ + +## Brute force + + +```octave +setdefaults +epsilon = 0.039; % kcal/mol +sigma = 2.934; % Angstrom +x=linspace(2.8,6,200); % bond length in Angstrom + +Ex = lennard_jones(x,sigma,epsilon); + +[Emin,imin]=min(Ex); + +plot(x,Ex,x(imin),Emin,'o') +``` + + +![svg](lecture_08_files/lecture_08_3_0.svg) + + + +```octave +x(imin-1) +x(imin) +x(imin+1) +``` + + ans = 3.2824 + ans = 3.2985 + ans = 3.3146 + + +## Golden Search Algorithm + +We can't just look for a sign change for the problem (unless we can take a derivative) so we need a new approach to determine whether we have a maximum between the two bounds. + +Rather than using the midpoint of initial bounds, here the problem is more difficult. We need to compare the values of 4 function evaluations. The golden search uses the golden ratio to determine two interior points. + +![golden ratio](goldenratio.png) + +Start with bounds of 2.5 and 6 Angstrom. + + +```octave +% define Au atomic potential +epsilon = 0.039; % kcal/mol +sigma = 2.934; % Angstrom +Au_x= @(x) lennard_jones(x,sigma,epsilon); + +% calculate golden ratio +phi = 1/2+sqrt(5)/2; +% set initial limits +x_l=2.8; +x_u=6; + +% Iteration #1 +d=(phi-1)*(x_u-x_l); + +x1=x_l+d; % define point 1 +x2=x_u-d; % define point 2 + + +% evaluate Au_x(x1) and Au_x(x2) + +f1=Au_x(x1); +f2=Au_x(x2); +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') +hold on; + +if f2x2 + plot(x4,f4,'*',[x1,x2],[f1,f2]) + x1=x2; + f1=f2; + else + plot(x4,f4,'*',[x3,x2],[f3,f2]) + x3=x2; + f3=f2; + end + x2=x4; f2=f4; +else + error('no minimum in bracket') +end +hold off +``` + + +![svg](lecture_08_files/lecture_08_11_0.svg) + + + +```octave +p = polyfit([x1,x2,x3],[f1,f2,f3],2); +x_fit = linspace(x1,x3,20); +y_fit = polyval(p,x_fit); + +plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o') +hold on +if f2x2 + plot(x4,f4,'*',[x1,x2],[f1,f2]) + x1=x2; + f1=f2; + else + plot(x4,f4,'*',[x3,x2],[f3,f2]) + x3=x2; + f3=f2; + end + x2=x4; f2=f4; +else + error('no minimum in bracket') +end +hold off + +``` + + +![svg](lecture_08_files/lecture_08_12_0.svg) + + + +```octave +p = polyfit([x1,x2,x3],[f1,f2,f3],2); +x_fit = linspace(x1,x3,20); +y_fit = polyval(p,x_fit); + +plot(x,Au_x(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o') +hold on +if f2x2 + plot(x4,f4,'*',[x1,x2],[f1,f2]) + x1=x2; + f1=f2; + else + plot(x4,f4,'*',[x3,x2],[f3,f2]) + x3=x2; + f3=f2; + end + x2=x4; f2=f4; +else + error('no minimum in bracket') +end +hold off + +``` + + +![svg](lecture_08_files/lecture_08_13_0.svg) + + +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) + +## Using the solutions to minimization for the nonlinear spring constant + +Now, we have two routines to find minimums of a univariate function (Golden Ratio and Parabolic). Let's use these to solve for the minimum energy given a range of applied forces to the single atom gold chain + +$E_{total}(\Delta x) = E_{LJ}(x_{min}+\Delta x) - F \cdot \Delta x$ + + +```octave +epsilon = 0.039; % kcal/mol +epsilon = epsilon*6.9477e-21; % J/atom +epsilon = epsilon*1e18; % aJ/J +% final units for epsilon are aJ + +sigma = 2.934; % Angstrom +sigma = sigma*0.10; % nm/Angstrom +x=linspace(2.8,6,200)*0.10; % bond length in um + +Ex = lennard_jones(x,sigma,epsilon); + +%[Emin,imin]=min(Ex); + +[xmin,Emin] = fminbnd(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6) + +plot(x,Ex,xmin,Emin,'o') +ylabel('Lennard Jones Potential (aJ/atom)') +xlabel('bond length (nm)') + +Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx; + +``` + + xmin = 0.32933 + Emin = -2.7096e-04 + + + +![svg](lecture_08_files/lecture_08_16_1.svg) + + + +```octave +N=50; +dx = zeros(1,N); % [in nm] +F_applied=linspace(0,0.0022,N); % [in nN] +for i=1:N + optmin=goldmin(@(dx) Etotal(dx,F_applied(i)),-0.001,0.035); + dx(i)=optmin; +end + +plot(dx,F_applied) +xlabel('dx (nm)') +ylabel('Force (nN)') +``` + + warning: Matlab-style short-circuit operation performed for operator | + warning: called from + goldmin at line 17 column 1 + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + warning: Matlab-style short-circuit operation performed for operator | + + + +![svg](lecture_08_files/lecture_08_17_1.svg) + + +For this function, it is possible to take a derivative and compare the analytical result: + + +```octave +dx_full=linspace(0,0.06,50); +F= @(dx) 4*epsilon*6*(sigma^6./(xmin+dx).^7-2*sigma^12./(xmin+dx).^13) +plot(dx_full,F(dx_full),dx,F_applied) +``` + + F = + + @(dx) 4 * epsilon * 6 * (sigma ^ 6 ./ (xmin + dx) .^ 7 - 2 * sigma ^ 12 ./ (xmin + dx) .^ 13) + + + + +![svg](lecture_08_files/lecture_08_19_1.svg) + + +## Curve-fitting +Another example is minimizing error in your approximation of a function. If you have data (now we have Force-displacement data) we can fit this to a function, such as: + +$F(x) = K_{1}\Delta x + \frac{1}{2} K_{2}(\Delta x)^{2}$ + + + +```octave +function SSE = sse_of_parabola(K,xdata,ydata) + % calculate the sum of squares error for a parabola given a function, func, and xdata and ydata + % output is SSE=sum of squares error + K1=K(1); + K2=K(2); + y_function = K1*xdata+1/2*K2*xdata.^2; + SSE = sum((ydata-y_function).^2); +end + +``` + + +```octave +[K,SSE_min]=fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]); +fprintf('\nThe nonlinear spring constants are K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\n',K) +fprintf('The mininum sum of squares error = %1.2e',SSE_min) +``` + + + The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.98 nN/nm^2 + The mininum sum of squares error = 7.35e-08 + + + +```octave +plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2) +``` + + +![svg](lecture_08_files/lecture_08_23_0.svg) + diff --git a/lecture_08/lecture_08_files/lecture_08_11_0.svg b/lecture_08/lecture_08_files/lecture_08_11_0.svg new file mode 100644 index 0000000..c61f790 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_11_0.svg @@ -0,0 +1,168 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.06 + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + gnuplot_plot_3a + + + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_12_0.svg b/lecture_08/lecture_08_files/lecture_08_12_0.svg new file mode 100644 index 0000000..b641032 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_12_0.svg @@ -0,0 +1,163 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + gnuplot_plot_3a + + + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_13_0.svg b/lecture_08/lecture_08_files/lecture_08_13_0.svg new file mode 100644 index 0000000..3e467e5 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_13_0.svg @@ -0,0 +1,163 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + gnuplot_plot_3a + + + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_16_1.svg b/lecture_08/lecture_08_files/lecture_08_16_1.svg new file mode 100644 index 0000000..ff27c29 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_16_1.svg @@ -0,0 +1,142 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.0004 + + + + + -0.0002 + + + + + 0 + + + + + 0.0002 + + + + + 0.0004 + + + + + 0.0006 + + + + + 0.2 + + + + + 0.3 + + + + + 0.4 + + + + + 0.5 + + + + + 0.6 + + + + + 0.7 + + + + + + + + + Lennard Jones Potential (aJ/atom) + + + + + bond length (nm) + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_17_1.svg b/lecture_08/lecture_08_files/lecture_08_17_1.svg new file mode 100644 index 0000000..1f3103b --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_17_1.svg @@ -0,0 +1,136 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0.0025 + + + + + -0.01 + + + + + 0 + + + + + 0.01 + + + + + 0.02 + + + + + 0.03 + + + + + 0.04 + + + + + + + + + Force (nN) + + + + + dx (nm) + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_19_1.svg b/lecture_08/lecture_08_files/lecture_08_19_1.svg new file mode 100644 index 0000000..2a891df --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_19_1.svg @@ -0,0 +1,135 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0.0025 + + + + + 0 + + + + + 0.01 + + + + + 0.02 + + + + + 0.03 + + + + + 0.04 + + + + + 0.05 + + + + + 0.06 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_23_0.svg b/lecture_08/lecture_08_files/lecture_08_23_0.svg new file mode 100644 index 0000000..8fa3bb1 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_23_0.svg @@ -0,0 +1,196 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + 0.0005 + + + + + 0.001 + + + + + 0.0015 + + + + + 0.002 + + + + + 0.0025 + + + + + 0 + + + + + 0.005 + + + + + 0.01 + + + + + 0.015 + + + + + 0.02 + + + + + 0.025 + + + + + 0.03 + + + + + 0.035 + + + + + 0.04 + + + + + + + + + gnuplot_plot_1a + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_3_0.svg b/lecture_08/lecture_08_files/lecture_08_3_0.svg new file mode 100644 index 0000000..2c38150 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_3_0.svg @@ -0,0 +1,147 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_6_1.svg b/lecture_08/lecture_08_files/lecture_08_6_1.svg new file mode 100644 index 0000000..ecf7a18 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_6_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_7_1.svg b/lecture_08/lecture_08_files/lecture_08_7_1.svg new file mode 100644 index 0000000..e63011c --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_7_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_8_1.svg b/lecture_08/lecture_08_files/lecture_08_8_1.svg new file mode 100644 index 0000000..873b04d --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_8_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lecture_08_files/lecture_08_9_1.svg b/lecture_08/lecture_08_files/lecture_08_9_1.svg new file mode 100644 index 0000000..28fd681 --- /dev/null +++ b/lecture_08/lecture_08_files/lecture_08_9_1.svg @@ -0,0 +1,169 @@ + + +Gnuplot +Produced by GNUPLOT 5.0 patchlevel 3 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -0.04 + + + + + -0.02 + + + + + 0 + + + + + 0.02 + + + + + 0.04 + + + + + 0.06 + + + + + 0.08 + + + + + 2.5 + + + + + 3 + + + + + 3.5 + + + + + 4 + + + + + 4.5 + + + + + 5 + + + + + 5.5 + + + + + 6 + + + + + + + + + gnuplot_plot_1a + + + + + + gnuplot_plot_2a + + + + + + gnuplot_plot_3a + + + + + + gnuplot_plot_4a + + + + + + gnuplot_plot_5a + + + + + + gnuplot_plot_6a + + + + + + + + + + + + + \ No newline at end of file diff --git a/lecture_08/lennard_jones.m b/lecture_08/lennard_jones.m new file mode 100644 index 0000000..d18a6ad --- /dev/null +++ b/lecture_08/lennard_jones.m @@ -0,0 +1,4 @@ +function E_LJ =lennard_jones(x,sigma,epsilon) + E_LJ = 4*epsilon*((sigma./x).^12-(sigma./x).^6); +end + diff --git a/lecture_08/octave-workspace b/lecture_08/octave-workspace new file mode 100644 index 0000000..990b83a Binary files /dev/null and b/lecture_08/octave-workspace differ