Biseccion.sci

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=
endfunction

// Estas variables controlan la convergencia
//   maxiter    número máximo de iteraciones
//   tol        tolerancia en el error relativo
maxiter=50;
tol=1e-3;

// definir booleano para indicar convergencia
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

// n es el contador de iteraciones
n=0;

while n<maxiter

    // Incrementar contador de iteraciones
    n=n+1;

    // Calcular xn y el valor de la función f(xn)
    xn=(a+b)/2;
    fxn=f(xn);

    // Checar el (poco probable) caso de que xn sea la solución
    if fxn==0 then
        e=0;
        converge=%T;
        break
    end

    // Identificar 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("o que 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 HUBO CONVERGENCIA")
    disp("Alcanzado máximo de "+string(maxiter)+" iteraciones")
end