%%
%-----------------------------------%
% Zeros of a function: use of fzero %
%-----------------------------------%
%-----------------------------%
% Ex. 3.0.9
clear all
close all
clc
% plot the function in order to be sure
% that there are the zeros where we are looking for
x=linspace(-5,5,200);
f=-3*x.^2+x.^4+x.^5/10;
g=0*x;
 
%plot the energy potential
hold on
%axis 'auto xy'
axis([-2.2  2  -3  2.5])
plot(x,f,'r-','LineWidth',2)
plot(x,g,'k')
legend('U(x)/E_c=-3x^2+x^4+x^5/10',...
    'Location','North')
xlabel('x','fontsize',16),ylabel('U(x)/E_c','fontsize',16)

%text(0,-0.8,'E','fontsize',16)
%text(0,-0.8,['\fontsize{16}\color{blue}E'])
hold off

%export the figure into an .eps file
%which later on can be inserted into a file
%(e.g., Word, LaTeX, ...)
%use command print
print -f -depsc 'Fig_poten.eps'
%N.B. if we want to print an additional Figure 2 use 
%print -f2 

%%
%Find root of continuous function of one variable with fzero
%first write an anonymous function
f = @(x)(-3*x^2+x^4+x^5/10);
%find the zero close to x=2
fzero(f,2)
%or find the zero in an interval [a,b]=[1,2]
xzero1=fzero(f,[1,2])

%similarly for a different interval
xzero2=fzero(f,[-3,-1])

%%
%---------------------------------------%
% Zeros of a function: bisection method %
%---------------------------------------%
%do the same with the `while' routine in 
%the interval [a,b]=[1,2]
%
clear all
a=1; fa=-3*a^2+a^4+a^5/10;
b=2; fb=-3*b^2+b^4+b^5/10;

while b-a>2*eps
    %define the middle point
    x=(a+b)/2;
    fx=-3*x^2+x^4+x^5/10;

    if (fa*fx)>0
        %eliminate left half interval
        a=x;
    else
        %eliminate right half interval
        b=x;
    end
end
xzero1=x

%%
%or equivalently
clear all
a=1; fa=-Inf;
b=2; fb=Inf;
while b-a>eps*b
   x=(a+b)/2;
   fx=-3*x^2+x^4+x^5/10;
  if sign(fx) == sign(fa)
    a=x; fa=fx;
  else
    b=x; fb=fx;
  end
end
xzero1=x
%%
%do the same with the `while' routine in 
%the interval [a,b]=[-3,-1]
%
clear all
a=-3; fa=-3*a^2+a^4+a^5/10;
b=-1; fb=-3*b^2+b^4+b^5/10;

while b-a>2*eps
    %define the middle point
    x=(a+b)/2;
    fx=-3*x^2+x^4+x^5/10;

    if (fa*fx)>0
        %eliminate left half interval
        a=x;
    else
        %eliminate right half interval
        b=x;
    end
end
xzero2=x

%%
%-----------------------------%
% Ex. 3.0.10
clear all
close all
clc
% plot the function in order to be sure
% that there are the zeros where we are looking for

x=linspace(0,5,200);
f=x.^3-7*x.^2+14*x-6;
g=0*x;
 
%plot the energy potential
hold on
%axis 'auto xy'
%axis([-2.2  2  -3  2.5])
plot(x,f,'r-','LineWidth',2)
plot(x,g,'k')
legend('f','Location','North')
xlabel('x','fontsize',16),ylabel('f','fontsize',16)

%text(0,-0.8,'E','fontsize',16)
%text(0,-0.8,['\fontsize{16}\color{blue}E'])
hold off

%%
a=0; fa=a^3-7*a^2+14*a-6;
b=1; fb=b^3-7*b^2+14*b-6;

while b-a>2*10^(-2)
    %define the middle point
    x=(a+b)/2;
    fx=x^3-7*x^2+14*x-6;

    if (fa*fx)>0
        %eliminate left half interval
        a=x;
    else
        %eliminate right half interval
        b=x;
    end
end
xzero=x

%compare with
f = @(x)x^3-7*x^2+14*x-6;
xzero=fzero(f,[0,1])

%%
%--------------------------------------------%
% Zeros of a function: Newton-Raphson method %
%--------------------------------------------%
%-----------------------------%
% Ex. 4.0.11
clear all
close all
clc
% Example of 1 iteration

a=1.5; b=3; N=100;
x=linspace(a,b,N);

f=exp(x)-x-4;
g=0*x;

% vector of derivative
%deltax=diff(x);
%deltaf=diff(f);
%fp=deltaf./deltax;

% first guess of the zero
val0=2.5;
% remembering that  x(i)=a+(i-1)*(b-a)/(N-1);
%i_val0=round( (val0-a)*(N-1)/(b-a) + 1 );
% check:
%x(i_val0)
% one way of evaluating the derivative could be to consider
%fp_val0=fp(i_val0)

% alternatively, and more simply, we can evaluate the derivative of f 
% in val0 only as
eps=10^(-10);
fp_val0=((exp(val0)-val0-4)-(exp(val0-eps)-(val0-eps)-4))/(eps)

% or even better exactly as
fp_val0_ex=exp(val0)-1;

f_val0=exp(val0)-val0-4;

%line of the first step of the Newton-Raphson method:
h=fp_val0*(x-val0)+f_val0;

y=linspace(0,f_val0,N);
val_vect=0*y+val0;

%second step point of the Newton-Raphson method:
val1=val0-f_val0/fp_val0

hold on
plot(x,f,'r-','LineWidth',2)
legend('f(x)=exp(x)-x-4',...
    'Location','North')
plot(x,g,'k')
plot(x,h,'b-')
plot(val_vect,y,'k--')
xlabel('x','fontsize',16),ylabel('f(x)','fontsize',16)
text(val0,-1,'x_0','fontsize',16)
text(val1,-1,'x_1','fontsize',16)
hold off

%export the figure into an .eps file
%which later on can be inserted into a file
%(e.g., Word, LaTeX, ...)
%use command print
print -f -depsc 'Fig_Newton-R-method.eps'
%print -f -depsc 'Fig_Newton-R-method.eps'
%%
%-----------------------------%
% Ex. 4.0.12
clear all
close all
clc
%Find root of continuous function of one variable with fzero
%first write an anonymous function
f = @(x)(exp(x)-x-4);
%find the zero close to x=2
fzero(f,2)

%%
%-----------------------------%
% Ex. 4.0.13
clear all
format long
% first guess
val0=2.5;
% counter
kount=0;
% initialise the loop while
f_val1=1;
while f_val1>2*eps
    f_val0=exp(val0)-val0-4;
    % derivative in val0
    fp_val0=exp(val0)-1;

    val1=val0-f_val0/fp_val0;
    f_val1=exp(val1)-val1-4;
    
    % the new guess now is val1
    val0=val1;
    kount=kount+1;
end
val0
f_val1
kount

%%
%-----------------------------%
% Ex. 4.0.14
clear all

% same function as previous excercise
a=1; fa=exp(a)-a-4;
b=2; fb=exp(b)-b-4;

kount=0;
while b-a>2*eps
    %define the middle point
    x=(a+b)/2;
    fx=exp(x)-x-4;
    if (fa*fx)>0
        %eliminate left half interval
        a=x;
    else
        %eliminate right half interval
        b=x;
    end
    kount=kount+1;
end
x
kount

%%
%------------------------------------------------------------%
% Newton-Raphson method to find the max or min of a function %
%------------------------------------------------------------%
%-----------------------------%
% Ex. 4.0.16
clear all
close all
a=-2; b=2; N=100;
x=linspace(a,b,N);
f=-3*x.^2+x.^4+0.1*x.^5;

deltax=diff(x);
deltaf=diff(f);

fp=deltaf./deltax;
il=(1:1:N-1); ir=(2:1:N);

g=0*x;

hold on
plot(x,f,'r-','LineWidth',2)
plot(x(il)+deltax/2,fp,'ks', x,-6*x+4*x.^3+x.^4/2,'b')
legend('f(x)=-3x^2+x^4+x^5/10','approximate derivative','f´(x)',...
    'Location','NorthWest')
plot(x,g,'k-')
xlabel('x','fontsize',16), ylabel('f´(x)','fontsize',16)
hold off
%print -f -depsc '/home/fmarchet/classes/computacion/mater_compl/unit3/derivative.eps'
%%
clear all
f_val1=1;

%first guess
val0=1;

kount=0;

while f_val1>2*eps
    
    % first derivative in val0
    f_val0=-6*val0+4*val0^3+val0^4/2;
    % second derivative in val0
    fp_val0=-6+12*val0^2+2*val0^3;

    val1=val0-f_val0/fp_val0;
    
    f_val1=-6*val1+4*val1^3+val1^4/2;
    
    % the new guess now is val1
    val0=val1;
    kount=kount+1;
end
val0
f_val1
kount

%%
clear all

a=1; fa=-6*a+4*a^3+a^4/2;
b=2; fb=-6*b+4*b^3+b^4/2;

kount=0;
while b-a>2*eps
    %define the middle point
    x=(a+b)/2;
    fx=-6*x+4*x^3+x^4/2;
    if (fa*fx)>0
        %eliminate left half interval
        a=x;
    else
        %eliminate right half interval
        b=x;
    end
    kount=kount+1;
end
x
kount
%%
%------------------------------------%
% Excercise 4.0.17                    %
%------------------------------------%
clear all
close all
clc

x0=2; y0=-3; R=5;
% define the two semicircles (note that they are meaningfull only for 
% R^2-(x-x0)^2 >= 0
% i.e., x<R+x0 and x>-R+x0
a=-R+x0
b=R+x0

c1=@(x)( y0+sqrt(R^2-(x-x0).^2) );
c2=@(x)( y0-sqrt(R^2-(x-x0).^2) );

f1=@(x)( exp(-x/2)-4-c1(x) )
f2=@(x)( exp(-x/2)-4-c2(x) )

x=linspace(a,b,100);
% let's plot to have an idea where the crossing points are
figure(1)
hold on
plot(x,f1(x),'-b')
plot(x,f2(x),'-r')
plot(x,0*x,'--k')
hold off

% first zero
a1=a; b1=0; kount1=0;
while abs(b1-a1)>2*eps
    %define the middle point
    xm1=(a1+b1)/2;
    if (f1(a1)*f1(xm1))>0
        %eliminate left half interval
        a1=xm1;
    else
        %eliminate right half interval
        b1=xm1;
    end
    kount1=kount1+1;
end
kount1
xsol1=xm1

% second zero
a2=5; b2=b; kount2=0;
while abs(b2-a2)>4*eps
    %define the middle point
    xm2=(a2+b2)/2;
    if (f2(a2)*f2(xm2))>0
        %eliminate left half interval
        a2=xm2;
    else
        %eliminate right half interval
        b2=xm2;
    end
    kount2=kount2+1;
end
kount2
xsol2=xm2

%%
% plotting again with final answer
% circle

th=linspace(0,2*pi,1000);
xc=x0+R*cos(th);
yc=y0+R*sin(th);
% exponential
x=linspace(-4,8,1000);
y=exp(-x/2)-4;

% values of crossing points
ysol1=exp(-xsol1/2)-4;
ysol2=exp(-xsol2/2)-4;

% plot
figure(2)
hold on
plot(xc,yc,'r','LineWidth',2)
plot(x,y,'g','LineWidth',2)
Legend('circle','exp(-x/2)-4')
plot(xsol1,ysol1,'ks','LineWidth',2,'MarkerFaceColor','b',...
                'MarkerSize',10)
plot(xsol2,ysol2,'ks','LineWidth',2,'MarkerFaceColor','b',...
                'MarkerSize',10)
xlabel('x','Fontsize',18),ylabel('y','Fontsize',18)
axis('equal')
hold off

print -f -depsc 'cross_circle_exp.eps'


%%
%---------------------%
% Energy conservation %
%---------------------%
%-----------------------------%
% Ex. 5.0.18
clear all

U=@(y)(-3*y.^2+y.^4+y.^5/10);
x=linspace(-3,3,100)
hold on
plot(x,U(x))
plot(x,0*ones(length(x),1),'k--')
axis([-3,3,-5,5])
% there are two roots (zeros)
% one between [1,2] and the other between [-3,-1]
a=1; b=2;
U(a)*U(b) % check it is negative
i=0; % contador
error=eps; %error=10^(-2);
while abs(b-a)>2*error
   xm=(a+b)/2;
   if U(a)*U(xm)>0
       a=xm;
   else
       b=xm;
   end
   i=i+1;
end
xm
U(xm)
i
plot(xm,U(xm),'rs')
hold off

%%
% exercise 5.0.18 with Newton-Raphson
clear all
close all
clc
U=@(y)(-3*y.^2+y.^4+y.^5/10);
dU=@(y)(-6*y+4*y.^3+y.^4/2);
% initial guess
x0=-7;
error=eps;
i=0;
while abs(U(x0))>2*error
   x1=x0-U(x0)/dU(x0); 
   x0=x1;
   i=i+1;
end
x1
i
x=linspace(-3,3,100);
hold on
plot(x,U(x))
plot(x1,U(x1),'rs')
plot(x,0*ones(length(x),1),'k--')
axis([-3,3,-5,5])
hold off