Biseccion.sce

clear
disp("Método de bisección para f(x)=0");

// Definir la función que se va a resolver
function y=f(x)
    y=cos(x)-x^2;
endfunction

// Estas variables controlan la convergencia
//   maxiter    número máximo de iteraciones
//   tol        tolerancia en el error relativo
//   converge   booleano para indicar si converge o no
maxiter=50;
tol=5e-8;
converge=%F;

disp("Intervalo de búsqueda [a,b]?");
a=input("a = ");
b=input("b = ");

// Calcular la función en ambos extremos
fa=f(a);
fb=f(b);

// Si ambas tienen el mismo signo, no se puede proceder
if (fa*fb)>0 then
    disp("f(a) y f(b) tienen el mismo signo");
    disp("No parece haber una raíz entre a y b");
    disp("  tal vez la funcion no sea continua o")
    disp("  el intervalo contiene más de una solución")
    abort
end

// Para el cálculo del error en la primera iteración, tomar
// como valor previo el de menor valor absoluto de f(x)
if abs(fa)<abs(fb) then
    xp=a;
else
    xp=b;
end

// Iteraciones (n es el contador)
n=0;
while n<maxiter
    
    n=n+1;
    
    // Calcular xn y el valor de la función f(xn)
    xn=(a+b)/2;
    fxn=f(xn);

    // Checar en cuál intervalo está la solución
    if (fa*fxn)<0 then
        // La solución está entre a y xn
        // Descartar b y remplazar por el valor de xn
        b=xn;
        fb=fxn;
    elseif (fb*fxn)<0 then
        // La solución está entre xn y b
        // Descartar a y remplazar por el valor de xn
        a=xn;
        fa=fxn;
    else
        // Cuando ninguno de los dos productos da negativo
        disp("No se pudo encontrar la solución,")
        disp("puede haber más de una solución")
        disp("que o la función no sea continua")
        abort
    end
    
    // Checar convergencia (error relativo aproximado)
    e=abs((xn-xp)/xn);
    if e<tol then
        // El error es menor que la tolerancia
        // Señalar que ya convergió y salir del ciclo
        converge=%T;
        break
    end

    // Guardar xn como x previa para siguiente iteración
    xp=xn;
end

// Reportar los resultados (si convergió)
if (converge) then
    disp("Converge en "+string(n)+" iteraciones")
    disp("Solución encontrada:")
    disp("    x = "+string(xn))
    disp(" f(x) = "+string(fxn))
    disp("error = "+string(e*100)+"%")
else
    disp("NO CONVERGE")
    disp("Alcanzado máximo de "+string(maxiter)+" iteraciones")
end