%% Week 8
% Matlab codes for week 8

%% Problem 1

format long
xstar=((27^(1/2)*5240^(1/2))/27 + 352/27)^(1/3) - 26/(9*((27^(1/2)*5240^(1/2))/27 + 352/27)^(1/3)) - 2/3

syms x
p=x^3+2*x^2+10*x-20
pder=diff(p)
p=matlabFunction(p);
pder=matlabFunction(pder);
fplot(p,[1,2])

%%
% b)
x=2, 
for i=1:6 
    i
    x=x-p(x)/pder(x) 
    disp('Press <enter>')
    pause
end

%%
% c)
x=2, for i=1:50 i, x=20/(x^2+2*x+10), disp('Press <enter>'), pause, end

%%
% d)

fzero(p,1,optimset('Display','iter'))
pause

fzero(p,[1,2],optimset('Display','iter'))

%% Problem 2

clear error
syms x
f=x^3-5*x^2+8*x-4;
fder=diff(f);
f=matlabFunction(f);
fder=matlabFunction(fder);

x=3, for i=1:15 i, x=x-f(x)/fder(x), error(i)=abs(x-2); disp('Press <enter>'), pause, end

plot(log(error(1:end-1)),log(error(2:end)))
axis equal

%% Problem 3

f=@(x)x.*(63*x.^4-70*x.^2+15)/8;
F=@(x)x-f(x)/10;

Fs=sym(F);
Fders=diff(Fs); Fder=matlabFunction(Fders);

fplot(F,[0.8,1]), axis equal, title('Graph of F'), hold on, fplot(@(x)x,[0.8,1])
pause, hold off
fplot(Fder,[0.8,1]), title('Graph of F''')
q=Fder(0.8)

x0=1, x1=F(x0)
k_est=log(5e-7/abs(x1-x0)*(1-q))/log(q)

pause

x=1, for k=1:112 k, x=F(x), end

%% Problem 4

xstar=2;
F=@(x)x-1+1/(x-1);
Fs=sym(F)
Fders=diff(Fs)
Fder=matlabFunction(Fders);
Fder(xstar) % This shows at least second order convergence

pause

fplot(F,[1.5,2.5]), hold on, fplot(@(x)x,[1.5,2.5]), title('Graph of F')
hold off

%% Problem 5

syms x y z 
Fs=[x^2+2*y^2-y-2*z;x^2-8*y^2+10*z;x^2-7*y*z]
Js=jacobian(Fs)
F=matlabFunction(Fs,'vars',{[x;y;z]})
J=matlabFunction(Js,'vars',{[x;y;z]})

x=[1;1;1],
for k=1:8 z=J(x)\F(x); x=x-z, end

pause

fsolve(F,[1;1;1],optimset('Display','iter'))

%% Problem 6

F=@(x)(1-x(1))^2+5*(x(2)-x(1)^2)^2
fminunc(F,[2,2])
