#Métodos Numéricos - 2023. IMERL, Fing
#Código de ejemplo para el ejercicio 3 del práctico 4
clear all;
clc;
pkg load symbolic;

#Elegir qué parte se quiere ejecutar
parte='c';

f=@(x) 816*x^3 + (-3835)*x^2 + 6000*x + (-3125);
f_der=@(x) 816*3*x^2-3835*2*x+6000;

###################################
#Parte a
###################################
#syms x;
#factor(816*x^3 + (-3835)*x^2 + 6000*x + (-3125));

#Las raíces:
r1=5/3;
r2=25/16;
r3=25/17;

###################################
#Parte b
###################################
if parte=='b'
  figure(1);
  X=linspace(1,3);
  Y=816 .*X.^3 + (-3835).*X.^2 + 6000 .*X + (-3125);
  plot(X,Y,'LineWidth',2,'-b');
  set(gca,'fontsize', 26);
  hold on;
  grid on;

  line ([1.43 1.71], [0 0], "linestyle", "--", "color", [1, 0, 1],'LineWidth',1.5);
  scatter([r1,r2,r3],zeros(3,1),'k','s');
  title('816x^3-3835x^2+6000x-3125','FontSize',28);
endif

###################################
#Parte c - Método de Newton
###################################
if parte=='c'
  x0=1.5;
  figure(2);
  x_new=x0*ones(1,1);
  kmax=20;
  k=1;

  figure(2);
  X=linspace(1.43,1.71);
  Y=816 .*X.^3 + (-3835).*X.^2 + 6000 .*X + (-3125);
  plot(X,Y,'LineWidth',2,'-b');
  grid on;
  hold on;
  set(gca,'fontsize', 26);
  scatter([r1,r2,r3],zeros(3,1),'k','s');
  title('Método de Newton','FontSize',28);
  punto=plot([x_new(k)],[0],"marker", "*", "color", [1, 0, 1],'markersize',10);

  pause(2);
  delete(punto);
  k=k+1;

  x_new=[x_new,x_new(k-1)-f(x_new(k-1))/f_der(x_new(k-1))];
  punto=plot([x_new(k)],[0],"marker", "*", "color", [1, 0, 1],'markersize',10);

  while (k<kmax) & (abs(x_new(end) - x_new(end-1)) > eps*abs(x_new(end)))
    pause(2);
    delete(punto);
    k=k+1;

    x_new=[x_new,x_new(k-1)-f(x_new(k-1))/f_der(x_new(k-1))];
    punto=plot([x_new(k)],[0],"marker", "*", "color", [1, 0, 1],'markersize',10);
  endwhile
endif

###################################
#Parte d - Método de la secante
###################################
if parte=='d'
  kmax=20;
  x=[1;2];
  k=3;

  figure(3);
  X=linspace(1.43,1.71);
  Y=816 .*X.^3 + (-3835).*X.^2 + 6000 .*X + (-3125);
  plot(X,Y,'LineWidth',2,'-b');
  set(gca,'fontsize', 26);
  hold on;
  grid on;
  line ([1.43 1.71], [0 0], "linestyle", "--", "color", [1, 0, 1],'LineWidth',1.5);
  scatter([r1,r2,r3],zeros(3,1),'k','s');
  title('Método de la secante','FontSize',28);
  linea_a=plot(x(end-1)*ones(100,1),linspace(-1,1,100),"linestyle", "--", "color", [1, 0,0],'LineWidth',2);
  linea_b=plot(x(end)*ones(100,1),linspace(-1,1,100),"linestyle", "--", "color", [0, 1, 0],'LineWidth',2);
  linea=line([x(k-1),x(k-2)],[f(x(k-1)),f(x(k-2))],'LineWidth',2);

  while (kmax>k) & (abs(x(end)-x(end-1)) > eps*abs(x(end)))
    pause(2);
    delete(linea);
    delete(linea_a);
    delete(linea_b);

    pend=(f(x(k-1))-f(x(k-2,1)))/(x(k-1,1)-x(k-2,1));
    x_=x(k-1)-f(x(k-1))/pend;
    x=[x;x_];
    k=k+1;

    linea=line([x(k-1),x(k-2)],[f(x(k-1)),f(x(k-2))],'LineWidth',2);
    linea_a=line([x(k-1),x(k-1)],[-3,2],'LineWidth',2,'linestyle','--','color',[0,0,0]);
    linea_b=line([x(k-2),x(k-2)],[-3,2],'LineWidth',2,'linestyle','--','color',[0,0,0]);
  end

endif

###################################
#Parte e - Método de bisección
###################################
if parte=='e'
  a=1;
  b=2;
  m=(a+b)/2;
  f_m = f(m);
  k=1;
  kmax=20;

  figure(4);
  X=linspace(1.43,1.71);
  Y=816 .*X.^3 + (-3835).*X.^2 + 6000 .*X + (-3125);
  plot(X,Y,'LineWidth',2,'-b');
  set(gca,'fontsize', 26);
  hold on;
  grid on;
  line ([1.43 1.71], [0 0], "linestyle", "--", "color", [1, 0, 1],'LineWidth',1.5);
  scatter([r1,r2,r3],zeros(3,1),'k','s');
  title('Método de bisección','FontSize',28);
  linea_a=plot(a*ones(100,1),linspace(-1,1,100),"linestyle", "--", "color", [1, 0,0],'LineWidth',2);
  linea_b=plot(b*ones(100,1),linspace(-1,1,100),"linestyle", "--", "color", [0, 1, 0],'LineWidth',2);

  while (kmax>k) & (b-a>eps)
    pause(2);
    delete(linea_a);
    delete(linea_b);

    f_a = f(a);
    f_b = f(b);

    if f_a*f_m>0
      a=m;
    else
      b=m;
    endif

    m=(a+b)/2;
    f_m=f(m);
    k=k+1;

    linea_a=plot(a*ones(100,1),linspace(-1,1,100),"linestyle", "--", "color", [1, 0,0],'LineWidth',2);
    linea_b=plot(b*ones(100,1),linspace(-1,1,100),"linestyle", "--", "color", [0, 1, 0],'LineWidth',2);

  endwhile


endif