updating with problem 4

README.md

@@ -152,6 +152,22 @@ print('figure01','-dpng')
 
 #Problem 3
 
+### Code
+
+``` MatLab
+
+%[root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)
+fun = @(x)(x-1)*exp(-(x-1)^2);
+d_fun = @(x) exp(-(x - 1)^2) - exp(-(x - 1)^2)*(2*x - 2)*(x - 1);
+root = zeros(1,5);
+ea = zeros(1,5);
+iter = zeros(1,5);
+for y = 1:5
+    [root(y),ea(y),iter(y)]=newtraph(fun,d_fun,1.2,.0001,y);
+end
+table = [iter' root' ea'];
+
+```
 ### Divergence of Newton-Raphson Method
 
 | iteration | x_i | approx error |
@@ -173,3 +189,73 @@ print('figure01','-dpng')
 | 3 | 1.0000 | 0.0011 |
 | 4 | 1.0000 | 0.0000 |
 | 5 | 1.0000 | 0.0000 |
+
+
+#Problem 4
+
+```MatLab
+epsilon = 0.039; % units are [kcal/mol]
+epsilon = epsilon*6.9477e-21; % [J/atom]
+epsilon = epsilon*1e18; % [aJ/J]
+% episilon ends up being in terms of aJ
+
+sigma = 2.934; % for Angstrom
+sigma = sigma*0.10; % nm/Angstrom
+
+%setting up LJ
+lennard_jones = @(x,sigma,epsilon) 4*epsilon*((sigma./x).^12-(sigma./x).^6);
+[x,E,ea,its] = goldmin(@(x) lennard_jones(x,sigma,epsilon),3.2,3.5)
+
+figure(1)
+parabolic(3.2,3.5)
+
+epsilon = 0.039; % [kcal/mol]
+epsilon = epsilon*6.9477e-21; % [J/atom]
+epsilon = epsilon*1e18; % [aJ/J]
+% epsilon ends up being in terms of aJ
+
+sigma = 2.934; % Angstrom
+sigma = sigma*0.10; % nm/Angstrom
+%finding bond length in [um]
+x=linspace(2.8,6,200)*0.10;
+ex = lennard_jones(x,sigma,epsilon);
+
+[xmin,emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)
+
+figure(2)
+plot(x,ex,xmin,emin,'o')
+ylabel('Lennard Jones Potential [aJ/Atom]')
+xlabel('Bond Length [nm]')
+title('LJ Potential vs Bond Length');
+
+e_total = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;
+
+N=30;
+warning('off')
+dx = zeros(1,N); % [in nm]
+F_applied=linspace(0,0.0022,N); % [in nN]
+for i=1:N
+    optmin=goldmin(@(dx) e_total(dx,F_applied(i)),-0.001,0.035);
+    dx(i)=optmin;
+end
+
+plot(dx,F_applied)
+xlabel('dx [nm]')
+ylabel('Force [nN]')
+title('Force vs dx')
+
+dx_full = linspace(0,0.06,N);
+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)
+
+[K,SSE_min] = fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]);
+%fprintf('\n Nonlinear spring constants = K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\n',K)
+%fprintf('The mininum sum of squares error = %1.2e \n',SSE_min)
+
+plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)
+```
+
+![Plot showing LJ vs Bond Length](./figures/figure02.png)
+![Plot showing Force vs dx](./figures/figure03.png)
+####Nonlinear spring constants = K1=%1.2f nN/nm and K2=%1.2f nN/nm^2
+####The mininum sum of squares error = %1.2e

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,3 @@
+function E_LJ =lennard_jones(x,sigma,epsilon)
+  E_LJ = 4*epsilon*((sigma./x).^12-(sigma./x).^6);
+end

parabolic.m

@@ -0,0 +1,41 @@
+function [E4,x4] = parabolic(x_l,x_u)
+epsilon = .039;
+sigma = 2.934;
+fun = @(x) 4*epsilon*((sigma./x).^12-(sigma./x).^6);
+x = 3.2933;
+x_l = 3.22;
+x_u = 3.45;
+
+x1 = x_l;
+x2 = mean([x_l,x_u]);
+x3 = x_u;
+
+f1 = fun(x1);
+f2 = fun(x2);
+f3 = fun(x3);
+p = polyfit([x1,x2,x3],[f1,f2,f3],2);
+x_fit = linspace(x1,x3,20);
+y_fit = polyval(p,x_fit);
+
+plot(x,fun(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')
+hold on
+if f2<f1 && f2<f3
+    x4=x2-0.5*((x2-x1)^2*(f2-f3)-(x2-x3)^2*(f2-f1))/((x2-x1)*(f2-f3)-(x2-x3)*(f2-f1));
+    f4=fun(x4);
+
+    if x4>x2
+        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')
+    f4 = 'no';
+end
+hold off
+E4=f4

prob3.m

@@ -0,0 +1,10 @@
+%[root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)
+fun = @(x)(x-1)*exp(-(x-1)^2);
+d_fun = @(x) exp(-(x - 1)^2) - exp(-(x - 1)^2)*(2*x - 2)*(x - 1);
+root = zeros(1,5);
+ea = zeros(1,5);
+iter = zeros(1,5);
+for y = 1:5
+    [root(y),ea(y),iter(y)]=newtraph(fun,d_fun,1.2,.0001,y);
+end
+table = [iter' root' ea'];

prob4.m

@@ -0,0 +1,59 @@
+epsilon = 0.039; % units are [kcal/mol]
+epsilon = epsilon*6.9477e-21; % [J/atom]
+epsilon = epsilon*1e18; % [aJ/J]
+% episilon ends up being in terms of aJ
+
+sigma = 2.934; % for Angstrom
+sigma = sigma*0.10; % nm/Angstrom
+
+%setting up LJ
+lennard_jones = @(x,sigma,epsilon) 4*epsilon*((sigma./x).^12-(sigma./x).^6);
+[x,E,ea,its] = goldmin(@(x) lennard_jones(x,sigma,epsilon),3.2,3.5)
+
+figure(1)
+parabolic(3.2,3.5)
+
+epsilon = 0.039; % [kcal/mol]
+epsilon = epsilon*6.9477e-21; % [J/atom]
+epsilon = epsilon*1e18; % [aJ/J]
+% epsilon ends up being in terms of aJ
+
+sigma = 2.934; % Angstrom
+sigma = sigma*0.10; % nm/Angstrom
+%finding bond length in [um]
+x=linspace(2.8,6,200)*0.10; 
+ex = lennard_jones(x,sigma,epsilon);
+
+[xmin,emin] = goldmin(@(x) lennard_jones(x,sigma,epsilon),0.28,0.6)
+
+figure(2)
+plot(x,ex,xmin,emin,'o')
+ylabel('Lennard Jones Potential [aJ/Atom]')
+xlabel('Bond Length [nm]')
+title('LJ Potential vs Bond Length');
+
+e_total = @(dx,F) lennard_jones(xmin+dx,sigma,epsilon)-F.*dx;
+
+N=30;
+warning('off')
+dx = zeros(1,N); % [in nm]
+F_applied=linspace(0,0.0022,N); % [in nN]
+for i=1:N
+    optmin=goldmin(@(dx) e_total(dx,F_applied(i)),-0.001,0.035);
+    dx(i)=optmin;
+end
+
+plot(dx,F_applied)
+xlabel('dx [nm]')
+ylabel('Force [nN]')
+title('Force vs dx')
+
+dx_full = linspace(0,0.06,N);
+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)
+
+[K,SSE_min] = fminsearch(@(K) sse_of_parabola(K,dx,F_applied),[1,1]);
+%fprintf('\n Nonlinear spring constants = K1=%1.2f nN/nm and K2=%1.2f nN/nm^2\n',K)
+%fprintf('The mininum sum of squares error = %1.2e \n',SSE_min)
+
+plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)