02_roots_and_optimization

Answers to Homework 2

Problem 2

part a

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('figure_1','-dpng')

part b

### Method Analysis Table
| solver | es | ea | iterations |
| --- | --- | --- | --- |
| falsepos | 100,1000 | 9.637e-06 | 202 |
| mod_secant | 100,1000 | 5.9066e-06 | 24 |
| bisect | 100,1000 | 4.1212e-06 | 8 |

part c

Problem 3

part a

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
table01 = [iter' root' ea'];

for y = 1:5
    [root(y),ea(y),iter(y)]=newtraph(fun,dfun,1.2,.0001,y);
end
table02 = [iter' root' ea'];

part b

### 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 |

part c

### 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),2.8,3.5)

figure(1)
parabolic(2.8,3.5)
print('figure_2','-dpng')

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 minimum sum of squares error = %1.2e \n',SSE_min)

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

Output:The nonlinear spring constants are K1=0.16 nN/nm and K2=-5.86 nN/nm^2 The minimum sum of squares error = 4.95e-08 x = 2.8000 E =-1.4347e-09 ea = 9.2100e-05 its = 24 E4 = -0.0390 ans =-0.0390 xmin = 0.3293 Emin = -2.7096e-04