02_roots_and_optimization

ME3255 HW2

#Problem 2

###Part A

function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
% bisect: root location zeroes
% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
% uses bisection method to find the root of func
% input:
% func = 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 func
% output:
% root = real root
% fx = function value at root
% ea = approximate relative error (%)
% iter = number of iterations
if nargin<3,error('at least 3 input arguments required'),end
test = func(xl,varargin{:})*func(xu,varargin{:});
if test>0,error('no sign change'),end
if nargin<4||isempty(es), es=0.0001;end
if nargin<5||isempty(maxit), maxit=50;end
iter = 0; xr = xl; ea = 100;
while (1)
  xrold = xr;
  xr = (xl + xu)/2;
  iter = iter + 1;
  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
  test = func(xl,varargin{:})*func(xr,varargin{:});
  if test < 0
    xu = xr;
  elseif test > 0
    xl = xr;
  else
    ea = 0;
  end
  if ea <= es || iter >= maxit,break,end
end
root = xr; fx = func(xr, varargin{:});

function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)
% bisect: root location zeroes
% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
% uses bisection method to find the root of func
% input:
% func = 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 func
% output:
% root = real root
% fx = function value at root
% ea = approximate relative error (%)
% iter = number of iterations
if nargin<3,error('at least 3 input arguments required'),end
test = func(xl,varargin{:})*func(xu,varargin{:});
if test>0,error('no sign change'),end
if nargin<4|isempty(es), es=0.0001;end
if nargin<5|isempty(maxit), maxit=50;end
iter = 0; xr = xl; ea = 100;
while (1)
  xrold = xr;
  xr = (xl + xu)/2;
  % xr = (xl + xu)/2; % bisect method
  xr=xu - (func(xu)*(xl-xu))/(func(xl)-func(xu)); % false position method
  iter = iter + 1;
  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
  test = func(xl,varargin{:})*func(xr,varargin{:});
  if test < 0
    xu = xr;
  elseif test > 0
    xl = xr;
  else
    ea = 0;
  end
  if ea <= es | iter >= maxit,break,end
end
root = xr; fx = func(xr, varargin{:});

function [root,ea,iter]=mod_secant(func,dx,xr,es,maxit,varargin)
% newtraph: Modified secant root location zeroes
% [root,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):
%   uses modified secant method to find the root of func
% input:
%   func = name of function
%   dx = perturbation fraction
%   xr = initial guess
%   es = desired relative error (default = 0.0001%)
%   maxit = maximum allowable iterations (default = 50)
%   p1,p2,... = additional parameters used by function
% output:
%   root = real root
%   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
iter = 0;
while (1)
  xrold = xr;
  dfunc=(func(xr+dx)-func(xr))./dx;
  xr = xr - func(xr)/dfunc;
  iter = iter + 1;
  if xr ~= 0
    ea = abs((xr - xrold)/xr) * 100;
  end
  if ea <= es || iter >= maxit, break, end
end
root = xr;

###Part B

solver	initial guess(es)	ea	number of iterations
falsepos	100,1000	9.6376e-06	202
mod_secant	100,1000	5.9066e-06	24
bisect	100,1000	4.1212e-06	8

###Part C

cat_cable = @(T) T/10.*cosh(10./T*30)+30-T/10-35;
[root,fx,ea,iter] = falsepos(cat_cable,100,1000,0.00001,10000);
[root1,fx1,ea1,iter1] = bisect(cat_cable,100,1000,0.00001,10000);
[root2,ea2,iter2] = mod_secant(cat_cable,100,1000,0.00001,10000);

%define T
T = root2;

%plotting the shape of the powerline
x = -10:0.1:50;
y = T/10.*cosh(10./T*x)+30-T/10;
%setDefaults
plot(x,y)
title('Final Powerline Shape')
xlabel('distance (m)')
ylabel('height (m)')
print('figure01','-dpng')

#Problem 3

Code

%[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
0	3	n/a
1	3.2857	8.6957
2	3.5276	6.8573
3	3.7422	5.7348
4	3.9375	4.9605
5	4.1182	4.3873

Convergence of Newton-Raphson Method

iteration	x_i	approx error
0	1.2	n/a
1	0.9826	22.1239
2	1.0000	1.7402
3	1.0000	0.0011
4	1.0000	0.0000
5	1.0000	0.0000

#Problem 4

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]);

plot(dx,F_applied,'o',dx,K(1)*dx+1/2*K(2)*dx.^2)