diff --git a/README.md b/README.md index 8e06c99..63c1bce 100644 --- a/README.md +++ b/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 diff --git a/figures/figure02.png b/figures/figure02.png new file mode 100644 index 0000000..c440ba6 Binary files /dev/null and b/figures/figure02.png differ diff --git a/figures/figure03.png b/figures/figure03.png new file mode 100644 index 0000000..9cc647f Binary files /dev/null and b/figures/figure03.png differ diff --git a/goldmin.m b/goldmin.m new file mode 100644 index 0000000..3dd1dd1 --- /dev/null +++ b/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/lennard_jones.m b/lennard_jones.m new file mode 100644 index 0000000..773b4b7 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/parabolic.m b/parabolic.m new file mode 100644 index 0000000..9499563 --- /dev/null +++ b/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 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') + f4 = 'no'; +end +hold off +E4=f4 \ No newline at end of file diff --git a/prob3.m b/prob3.m new file mode 100644 index 0000000..10703ad --- /dev/null +++ b/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']; \ No newline at end of file diff --git a/prob4.m b/prob4.m new file mode 100644 index 0000000..f91681c --- /dev/null +++ b/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) \ No newline at end of file