% ==========
% Apartado I
% ==========
%
% Definición de funciones simples
%
f=inline('sqrt(exp(x))')
p=f(4)  % Evaluacion de la funcion f en el punto 4
X=linspace(-1,4);Y=f(X);plot(X,Y) % Representacion en [-1,4]
% ============
% Apartado III
% ============
%
% Definición de funciones simbolicas
%
clear all
syms x real
f=sqrt(exp(x))
p=subs(f,x,4)    % Evaluacion de la funcion f en el punto 4
ezplot(f,[-1,4]) % Representacion en [-1,4]
clear f
% ===========
% Apartado IV
% ===========
%
% Resolucion utilizando solve
%
% * * * EJERCICIO 1 * * *
f1=sin(x)+0.8*cos(x)
figure(1),ezplot(f1,[2,3]),grid
r1=solve(f1) % Solucion +k*pi
r1=solve('sin(x)+0.8*cos(x)') % Solucion +k*pi
figure(2),ezplot(f1),grid,pause
%
f2=x.^2-4*x+3.5-log(x)
figure(1),ezplot(f2,[1,3]),grid
r2=solve(f2) % Raices en [1,2] y [3,4]
figure(2),ezplot(f2),grid,pause
%
f3=(x-2.1).^2-7*x*cos(x)
figure(1),ezplot(f3,[1,2]),grid
r3=solve(f3) % Raices en [0,1],[1,2],[4,6],[6,8]
figure(2),ezplot(f3),grid,pause
% ==========
% Apartado V
% ==========
%
% Definición de T21_bisec
% Definicion de la función (debe estar en el fichero f1.m)
% ##################### Inicio del fichero ####################
function [c,err]=T21_bisec(f,a,b,tol,niter)
%
%   function [c,err]=T21_bisec(f,a,b,tol,niter)
%Resuelve una ecuacion no lineal utilizando biseccion
%
%Entrada
%	f	    nombre de la funcion (pasado como caddena de caracteres-fichero)
%	a,b	    extremos del intervalo
%	tol     tolerancia de la raiz y de la función respectivamente
%   niter   numero de iteraciones
%Salida
%	c		raiz aproximada
%   err	    error estimado para la raiz y de la función respectivamente
%
% Ejemplo:  [sol,errores]=T21_bisec(F1,-3,-2,0.001) % F1 con inline

if (nargin<5),niter=1000;end                  % Verificacion de los argumentos de entrada
if (nargin<4),tol=[eps eps];end                 
if (length(tol)~=2),tol=[tol(1) tol(1)];end
if nargin<3,error('ERROR T21_bisec: Numero insuficiente de argumentos');end
%
% Valoración de la función en los extremos
ya=feval(f,a);        
yb=feval(f,b);

if (ya*yb>0),error('ERROR T21_bisec: la funcion no cambia de signo');end
niter=min(niter,max(1,ceil((log(b-a)-log(tol(1)))/log(2)))); % Calculo del numero maximo de iteraciones

for k=1:niter
    c=(a+b)/2;% Calculo del nuevo valor
    yc=feval(f,c);                  
    if abs(yc) < tol(2);break;end;          % Valor de la función menor que la tolerancia
    if (sign(yc)==sign(yb))          % Cambio de los extremos
        b=c;yb=yc;
    else
        a=c;ya=yc;
    end
    if abs(b-a)<tol(1),break,end
end
if k>niter,disp('Agotadas las iteraciones sin convergencia');k=niter;end
if nargout>1, err=[abs(b-a) yc];end    % Errores generados
% ##################### fin del fichero ####################
% Resolución de las funciones
%
% * * * EJERCICIO 2 * * *
F1=vectorize(f1)
[c1b,err1b]=T21_bisec(@F1,2,3,1E-6)
F2=vectorize(f2)
[c2b,err2b]=T21_bisec(@F2,1,3,1E-6)
F3=vectorize(f3)
[c3b,err3b]=T21_bisec(@F3,1,2,1E-6)
%
% * * * EJERCICIO 3 * * *
% Modificación para obtener las iteraciones
% Cambiar la línea de definición de la función por
function [c,err,k]=T21_bisec(f,a,b,tol,niter)
%
% Modificación para obtener las resultados intermedios
% Insertar despues de la linea que calcula el nuevo valor de yc
disp([a,ya,b,yb,c, yc])
%
% * * * EJERCICIO 4 * * *
% Modificación de Réguila Falsi
% Cambiar la linea que calcula el nuevo valor c por 
c=a-ya*(b-a)/(yb-ya); % Calculo del nuevo valor
% Y eliminar la linea que calcula maxiter
%
% Resolucion utilizando T21_rfalsi
%
[c1r,err1r]=T21_rfalsi(F1,2,3,1E-6,100)
[c2r,err2r]=T21_rfalsi(F2,1,3,1E-6,100)
[c3r,err3r]=T21_rfalsi(F3,1,2,1E-6,100)
%
% Resolucion utilizando fzero
%
f=inline('x.^3-1');
options=optimset('Display','iter','TolX',1e-5,'MaxIter',100)
raiz=fzero(f,2,options)
%
% * * * EJERCICIO 5 * * *
options=optimset('Display','iter','MaxFunEvals',100,'TolFun',1e-10,'TolX',1e-5,'MaxIter',100)
R1=fzero(@F1,2,options)
% tambien resuelve fzero(f,2,options) siendo
%  f='sin(x)+0.8*cos(x)'
%  f=inline('sin(x)+0.8*cos(x)')
%  f=vectorize(f1)  % ES LA MEJOR OPCION
%
R2=fzero(F2,2,options)
R3=fzero(F3,2,options)
options=optimset
% ============
% Apartado VII
% ============
%
% Ecuaciones polinomicas
%
clear all
p=[1,2,3,4] % Define P(x)=x3+2x2+3x+4
roots(p)
P=poly2sym(p)
solve(P)
%
% * * * EJERCICIO 5 * * *
p=[1,0,2,-1,-3]
roots(p)
pause
% Representación del hundimiento de una esfera para diferentes densidades
p=[1,-3,0,0];
H=(03:0.1:0.7);N=length(H); % Para la resulución según el guion
%H=linspace(0,1);N=length(H); % Parar una representación mejor
R=zeros(size(H));
for i=1:N
    p(4)=4*H(i);
    r=roots(p)
    k=find(r>=0&r<=2)
    R(i)=max(r(k));
end
plot(H,50*R);grid on
xlabel('densidad');ylabel('porcentaje de diámetro hundido');
title('Hundimiento vs densidad en una esfera');