completed assignment

README.md

@@ -1,3 +1,61 @@
 # Roots and Optimization assignment
 
 Demo of uploading assignment
+
+## Solution to Gold Chain nonlinear spring constants
+
+!(Gold chain TEM image)[au_chain.jpg]
+
+!(Gold chain model (artist rendition))[Auchain_model.png]
+
+# Homework #3 plots of Gold chain F vs dx and Lennard-Jones Potential
+
+```matlab
+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)
+
+h1=figure(1)
+plot(x,Ex,xmin,Emin,'o')
+ylabel('Lennard Jones Potential (aJ/atom)')
+xlabel('bond length (nm)')
+saveas(h1,'potential_energy.png')
+
+Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;
+
+% Now with xmin determined find F vs dx for gold chain model
+% with 50 steps from 0 to 0.0022 nN
+
+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
+
+h2=figure(2)
+plot(dx,F_applied)
+xlabel('dx (nm)')
+ylabel('Force (nN)')
+saveas(h2,'force_vs_dx.png')
+```
+
+This script in included as `gold_chain_script.m`
+
+Output is two plots, 'potential_energy.png' and 'force_vs_dx.png'
+
+!(Lennard-Jones potential energy for no applied force)[potential_energy.png]
+
+!(Force vs displacement for gold chain)[force_vs_dx.png]

gold_chain_script.m

@@ -0,0 +1,40 @@
+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)
+
+h1=figure(1)
+plot(x,Ex,xmin,Emin,'o')
+ylabel('Lennard Jones Potential (aJ/atom)')
+xlabel('bond length (nm)')
+saveas(h1,'potential_energy.png')
+
+Etotal = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;
+
+% Now with xmin determined find F vs dx for gold chain model
+% with 50 steps from 0 to 0.0022 nN
+
+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
+
+h2=figure(2)
+plot(dx,F_applied)
+xlabel('dx (nm)')
+ylabel('Force (nN)')
+saveas(h2,'force_vs_dx.png')
+

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{:});

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
+

setdefaults.m

@@ -0,0 +1,3 @@
+set(0, 'defaultAxesFontSize', 16)
+set(0,'defaultTextFontSize',14)
+set(0,'defaultLineLineWidth',3)