02_roots_and_optimization

#Problem 2

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;

cat_cable= @(T) T/10.*cosh(10./T*30)+30-T/10-35;
[root,fx,ea,iter]=falsepos(cat_cable,100,1000,.00001,10000);
[root1,fx1,ea1,iter1]=bisect(cat_cable,100,1000,.00001,10000);
[root2,ea2,iter2]=mod_secant(cat_cable,100,1000,.00001,10000);
T=root2;
x=-10:.1:50;
y=T/10.*cosh(10./T*x)+30-T/10;
setDefaults
plot(x,y)
title('Final Shape of Powerline')
xlabel('Distance in meters')
ylabel('Height in meters')
print('figure01','-dpng')

###Method Analysis Table

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

#Problem 3

fun = @(x)(x-1)*exp(-(x-1)^2);
dfun = @(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,dfun,3,.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; % 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

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

[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)')

Etotal = @(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) Etotal(dx,F_applied(i)),-0.001,0.035);
    dx(i)=optmin;
end

plot(dx,F_applied)
xlabel('dx (nm)')
ylabel('Force (nN)')

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('\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 \n',SSE_min)

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

The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.86 nN/nm^2 The mininum sum of squares error = 4.95e-08