Extracted and adapted from http://www.ugr.es/~pmrivas/calculo-grado-en-informatic/bucles-biseccion-2014.wxmx

bucles-biseccion-2014.wxm

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 22.04.0 ] */
/* [wxMaxima: title   start ]
TEMA 5: Métodos numéricos de 
resolución de ecuaciones
   [wxMaxima: title   end   ] */


/* [wxMaxima: section start ]
Introducción al análisis numérico
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Hablamos de números máquina cuando nos referimos a los
números que tiene guardados el programa. Si introducimos un nº exacto 
wxMaxima elige el que más se aproxima a él. Con lo que está "redondeando".
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Errores de redondeo:
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
A lo anterior hay que añadir que cuando se pasa de sist. decimal a binario y se vuelve
al decimal, se acumulan los errores de redondeo y pueden pasar cosas como estas:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
fpprec:20;
bfloat(0.1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Sin embargo:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
bfloat(1/10);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Aritmética de ordenador
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
 Sabemos que el ordenador puede trabajar con números muy grandes
o muy pequeños; pero, si trabaja con ambos simultáneamente, el 
número pequeño puede perder su valor y hacerse cero debido al error 
de redondeo. Por eso hay que tener cuidado y recordar qué 
propiedades usuales en la artimética real (asociatividad, elemento neutro)
no son ciertas en la aritmética de ordenador.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Elemento neutro
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Tomamos un número muy pequeño, pero distinto de cero y vamos a ver cómo maxima 
interpreta que es cero:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
h:2.22045*10^(-17);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is(1.0+h=1.0);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Nos responde que es cierto que h+1=1, luego h sería cero.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Por encima, con números muy grandes puede hacer cosas raras.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
g:15.0+10^(20);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is(g-10^(20)=0);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
g-10^(20);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Aqui no sale igual, pero si los restáis cree que la diferencia es cero.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Cálculo de dos expresiones matemáticamente equivalentes
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
A la hora de realizar cálculos, la pérdida de precisión puede afectar al resultado.
Veamos algún ejemplo.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
remvalue(all);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Fijamos precisión 40 y consideramos el número a:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
fpprec : 40;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
a:bfloat(1-(10)^(-30));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Ahora vamos a calcular: a+1    y    (a^2-1)/(a-1). 
Deberían ser iguales:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
b:a+1;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
a^2-1;
a-1;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
c:(a^2-1)/(a-1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is(b=c);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Este es el resultado del efecto de cancelación de cifras significativas que tiene
lugar cuando se restan dos cantidades muy pequeñas.
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
Resolución numérica de ecuaciones con Maxima
   [wxMaxima: section end   ] */


/* [wxMaxima: subsect start ]
Comandos de resolución aproximada de ecuaciones
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Recordar allroots, realroots, bfallroots, nroots, find_root (Calcular raíz)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
allroots(x^5-x^3+6*x-81);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
nroots(x^5-x^3+6*x-81,-3,3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Comentamos ahora el comando: Ecuaciones --> Calcular raíz.
Se trata del comando find_root (que ya vimos en sesiones anteriores)
 es una aplicación del teorema de los 
ceros de Bolzano. 
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Por ejemplo, vamos a resolver la ecuación e^x+log(x)=0
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
f(x):=exp(x)+log(x);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
solve(f(x)=0,x);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
allroots(f(x));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Hemos visto cómo los comandos usuales de resolución de ecuaciones
 no funcionan.
Dibujamos la gráfica de f:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
wxplot2d([f(x)], [x,0,2],[y,-2,6])$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Parece que hay un cero entre 0 y 1.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
find_root(f(x), x, 0,1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Pero Maxima no nos da respuesta ya que ha intentado evaluar en cero
la función f. Hay que destacar que Maxima verifica en primer lugar 
si se dan las codiciones del T. de Bolzano; es decir, evalúa la función 
en los extremos que se le propones y comprueba que hay cambio de signo:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
find_root(f(x), x, 1,2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Ahora ya sí:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
find_root(f(x),x,0.1,1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Recomendación: Hacer los ejercicios: 
5.1, 5.2 y 5.3 de 
la página 88.
   [wxMaxima: subsect end   ] */


/* [wxMaxima: section start ]
Hágalo usted mismo: Breves conceptos 
de 
programación
   [wxMaxima: section end   ] */


/* [wxMaxima: subsect start ]
Bucles for
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Un bucle es un proceso repetitivo que se realiza un cierto número de veces.
Por ejemplo: Vamos a calcular los primeros 10 múltiplos de 7
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Ejemplo 1: Calcular los 10 primeros múltiplos de 7.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
for i:1  thru 10 do print(7*i)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Ejemplo 2: Calcular la suma de los 100 
primeros cuadrados.
Hay que añadir una variable adicional 
para guardar las sumas 
que vaya calculando.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
suma:0$
for i:1  thru 100 do suma:suma+i^2$
print(suma)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
if numer#false then numer:false else numer:true;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Ejemplo 3: Sumemos la serie geométrica de razón 
x=1/2
¡¡Cambiar a numérico!!!
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
numer:true;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
suma:1;
for i:1 thru 10 do suma:suma+1/(2^i)$
print(suma);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Si sumamos hasta 100 sumandos:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
suma:1.0;
for i:1 thru 100 do suma:suma+1/(2^i);
print(suma);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Hay un comando en Maxima para hacer sumatorias: sum
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
sum(1/(2^k), k, 1, 10), simpsum;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Bucles while, bucles unless, bucles do:
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
while condicion do expr               bucle while
unless condicion do expr              bucle unless
do expr                               bucle for
return (var)                          bucle for
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Utilizamos este tipo de bucles while cuando no sepamos cuántos pasos
hay que hacer, pero sí sepamos la condición que queremos que 
se cumpla en cada paso (o bucles unless, que se ejecuta siempre que
no se de una determinada condición). Así el bucle while acabará en cuanto dicha 
condición deje de darse (o en el bucle unless, en cuanto se dé dicha 
condición)
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Ejemplo 4: Queremos calcular cos(x)
comenzando en x = 0 e ir aumentando de 0.3 en 0.3 
hasta que el coseno deje de ser positivo.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
a:0$
while cos(a)>0 do (
                  print([a,cos(a)]),
                  a:a+0.3)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Y efectivamente, el siguiente es negativo:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
cos(1.8);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Condicionales
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
La segunda sentencia es la orden condicional if.
 Esta sentencia comprueba si se 
verifica una condición. 
Después, si la condición es verdadera, 
Maxima ejecutará una expresión1, 
y si es falsa ejecutará otra expresión2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
for i:1 thru 10 do(
   if log(i)<2 then 
   print("el logaritmo de ", i , 
   " es menor que 2") 
   else return(j:i)
                   )$
print("el logaritmo de ", j ,
" es mayor que 2")$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Observación: La variable contador 
i es local al bucle, así
que si queremos rescatarla en el return, hay que guardarla en 
una nueva variable j. 
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
EJERCICIOS A ENTREGAR el 
17 Diciembre: 
5.4 y 5.6 y UNO de estos tres:
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Ejercicio Propuesto 1: Escribir con un bucle while desde 1, 
y en intervalos de 0.5 en 0.5, 
los valores del logaritmo, siempre que éste sea menor que 2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Ejercicio Propuesto 2: Empezando en 7, y decreciendo con step 0.5, 
escribe [i,log(i)] siempre que log(i)>=1, con un bucle 
unless.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Ejercicio Propuesto 3: Utilizar la estructura if-then-else en 
un bucle for para,
empezando en 10, y decreciendo de uno en uno
hasta el 1, dar de salida las exponenciales que sean mayores 
que 1000 y parar cuando sea menor que 1000.
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
Método de Bisección
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
El método de bisección es una de las formas más elementales de 
buscar una solución de una ecuación. 
Se trata de una aplicación del teorema de Bolzano.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Un ejemplo particular
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
En primer lugar, presentamos la función
a la que le vamos a calcular el cero,
los extremos del intervalo
donde vamos a localizarlo y verificamos
la condición del teroema de Bolzano:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
f(x):=x^6+x-5;
a:0.0;
b:2.0;
f(a)*f(b);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Dibujamos la gráfica de f:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
wxplot2d([f(x)], [x,0,2],[y,-1,4]);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
¡Observad cómo ya de entrada trabajamos 
en numérico!
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Hacemos un bucle para localizar el cero de f,
dividiendo el intervalo [a,b] por la mitad, y
 lo hacemos en 10 pasos.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
for i:1 thru 10 do(
     c:(a+b)/2,
     if f(a)*f(c)<0 then b:c else a:c,
     print([a,b])
                  )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Ya sabemos que el cero tiene que estar 
"encerrado" en el último intervalo.
Podríamos seguir para asegurarnos del tercer 
decimal, pero mejor controlamos el error 
obligando a finalizar el bucle cuando 
el error sea menor que lo que hayamos fijado.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Control del error
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Si presentamos err como el error que prefijamos, 
el número de pasos que hay que hacer será: 
ceiling(log_2((b-a)/err)= el menor entero mayor 
o igual que log_2((b-a)/err)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
f(x):=x^6+x-5$
a:0.0$
b:2.0$
f(a)*f(b);
err:10^(-5)$
log2(x):=log(x)/log(2)$
pasos:ceiling(log2((b-a)/err));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
for i:1 thru pasos do(
     c:(a+b)/2,
if f(a)*f(c)<0 then b:c else a:c,
     print([a,b])
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
¿Y si hay suerte?
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Podría ocurrir que en un determinado paso, 
antes de llegar al final
del bucle, encontráramos el cero de f. 
Entonces, tendríamos que 
establecer cómo acabar si ocurriera esto.
Pero no podemos ir preguntando en cada paso 
is(f(c)=0) puesto que 
en numérico
esta pregunta puede tener respuestas sorprendentes.
Así que establecemos lo que sería un "cero gordo"; 
es decir, si llamamos prec a una precisión, 
por ejemplo, 10^(-5), en cuanto c  verifique
que |f(c)|<prec, eso será considerado como 
el cero ("gordo") y finalizará el bucle.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
f(x):=x^6+x-5$
a:0.0$
b:2.0$
f(a)*f(b)$
err:10^(-6)$
log2(x):=log(x)/log(2)$
prec:10^(-6)$
pasos:ceiling(log2((b-a)/err));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
for i:1 thru pasos do(
     c:(a+b)/2,
     if abs(f(c))<prec then (print("La solucion
     es exacta"),return(c))
     else 
     if f(a)*f(c)<0 then b:c else a:c)$
print("La solucion es ", c)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
kill(all);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Si queremos saber en qué paso se obtiene la solución:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
f(x):=x^6+x-5$
a:0.0$
b:2.0$
f(a)*f(b)$
err:10^(-6)$
log2(x):=log(x)/log(2)$
prec:10^(-5)$
pasos:ceiling(log2((b-a)/err));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
for i:1 thru pasos do(
     c:(a+b)/2,
     if abs(f(c))<prec then (print("La solucion
     es exacta"),j:i,return(c))
    else 
 if f(a)*f(c)<0 then b:c else a:c)$
print("La solucion es ", [c,j])$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Hemos obtenido la solución con la precisión y el error establecido en la salida 
i=17.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Fijaos que si aplicamos el comando find_root,
hasta el quinto decimal coinciden. 
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
find_root(f(x),x,a,b);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
¿SE PUEDE MEJORAR? Podéis seguir las indicaciones
de la página 59 del pdf.
   [wxMaxima: comment end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 22.04.0"$