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