La Bisección Con Python

por | 3 agosto, 2008

Hola, me he decidido este verano por aprender Python por varias razones:

  1. Con este lenguaje he podido encontrar una forma fácil de poder hacer que el usuario introduzca la expresión de una función matemática mediante el uso de eval, la única pega es que todavía no he encontrado una que calcule la derivada de la función introducida. Esto del eval se lo debo su conocimiento a Zootropo, el cual me atendió muy bien. ¡Gachiasss tioo!!
  2. He podido aprender a dibujar funciones matemáticas fácilmente con las librerías Matplotlib para Python, verdaderamente dan mucho juego y son muy fáciles de utilizar. Otra alternativa no tan espectacular son las VisualPython, pero con el mismo cometido, dar apoyo para dibujar gráficas matemáticas.
  3. Implementación de muchos métodos matemáticos con las librerías NumPy, como cálculo matricial, funciones matemáticas de origen, transformadas rápidas de Fourier, aplicación a la estadística y la aproximación numérica, etc..
  4. La unión de las tres anteriores razones hacen que obtengamos una especie de MatLab algo rústico bajo Python, pero el cual nos da una total libertad para investigar nuestros propios proyectos científicos pasando previamente por una programación bajo python, dándonos una casi total portabilidad de nuestros proyectos entre los diferentes tipos de sistemas operativos. Lo malo de Python en el cálculo matemático es que parece ser que no se lleva bien con grandes cálculos, aunque esto todavía tengo que experimentar en ello, recordad que estoy empezando a aprender Python todavía

Y bueno, para aprender Python me he decidido por echarme a la piscina olímpica y empezar con algo que no sea el típico Hola Mundo o los aburridos ejemplos con casos concretos de utilización de listas, diccionarios, bucles y demás; es que si no me aburro y opto por abandonar. A mi me mola generalizar y hacer algo con algo de chicha (( carne o jugo, con contenido )) que la vida ya es de por sí demasiado monótona Así que me he decidido por hacer un programa que calcule las raíces aproximadas de una función mediante el método de la bisección, dibujando la función matemática y una nube de puntos de los distintos valores aproximativos que se van obteniendo al aplicar la bisección.

El programa lo he dividido en tres ficheros con extensión .py, y que son:

  1. Metodos-> En el cual pedimos los datos necesarios por terminal y ejecutamos el método. Este es mi fichero main (principal) del proyecto. Lo he llamado Metodos porque al proyecto le voy a añadir en un futuro muy próximo los otros métodos: Secante, Regula-Falsi, Newton, Steffenssen y alguno más que vaya leyendo por ahí y me interese. Ya os informaré de ello conforme los vaya introduciendo.
  2. evaluar-> Este fichero va a ser muy importante en nuestro proyecto, ya que implementaremos un uso adecuado de eval y pediremos al usuario que nos introduzca la expresión analítica de la función matemática.
  3. Bisección-> En este fichero implementaremos dos funciones, una que nos hará la bisección y otra que nos dibujará la gráfica y la nube de puntos aproximativos del método. Ya veréis qué chulo nos queda!!!

Bueno vamos allá, el código de Metodos.py es el siguiente:

  1. #!/usr/bin/python
  2. # coding: latin-1
  3.  
  4. import os, sys
  5. import Biseccion
  6. from Numeric import *
  7.  
  8. a = float(input("Introduce el extremo inferior del intervalo\n"))
  9. b = float(input("Introduce el extremo superior del intervalo\n"))
  10. TOL = float(input("Introduce la tolerancia del metodo\n"))
  11. Biseccion.biseccion(a, b, TOL, 100)
  12. Biseccion.dibujar(a,b,TOL, 100)

Vamos a explicarlo un poquito, ¿vale?. Las tres primeras líneas son:

  1. #!/usr/bin/python
  2. # coding: latin-1
  3. import os, sys

Con la primera hacemos que bajo Lionux nuestro proyecto se pueda ejecutar dando doble clic al fichero principal, en nuestro caso Metodos.py, recordad que debéis antes darle permisos de ejecución. Las otras dos líneas sirven para poder introducir caracteres latinos como acentos y demás, como también que el programa se adecue a nuestro OS. Siempre es aconsejable poner estas líneas.
Luego tenemos dos líneas de importación que son:

  1. import Biseccion
  2. from Numeric import *

La primera nos va a permitir poder tener acceso a nuestro fichero Biseccion.py, así como poder acceder a sus funciones y clases implementadas en dicho fichero, mediante su correspondiente llamada, ahora luego os pongo el código de este fichero. Con la segunda línea vamos a importar la librería Numpy entera ( eso se debe al * ) si no queremos importarla en su totalidad podemos especificar mediante comas y después de la palabra clave import las cosas a importar de dicha librería, supongo que esto hará que el programa utilice menos recursos. Seguidamente tenemos unas líneas de petición de datos por terminal mediante las líneas:

  1. a = float(input("Introduce el extremo inferior del intervalo\n"))
  2. b = float(input("Introduce el extremo superior del intervalo\n"))
  3. TOL = float(input("Introduce la tolerancia del metodo\n"))

Aquí empezamos a ver algunas diferencias entre Python y, por ejemplo, C/C++, y es que la declaración explícita de variables con su especificación de tipo no es necesaria, esto es algo de lo que menos me agrada en Python, supongo que por ser matemático me gusta tener un cierto orden y rigor en algunas cosas, y esta es una de ellas; también creo que se debe a que en C++`se hizo un gran adelanto con el objeto cin en detrimento del scanf de C. En fin, me acostumbraré. Volviendo a estas tres últimas líneas lo que hacemos es pedir 3 datos por terminal al usuario final del programa, la a y la b son el extremo inferior y superior, respectivamente, del intervalo de defininición, y TOL es lo que los matemáticos llamamos tolerancia o cota de error que permitiremos al método de la bisección. Ambas tres, las pedimos mediante input y las transformamos debidamente a float, es decir; a reales en lenguaje matemático.
Terminamos el código del fichero con estas dos llamadas:

  1. Biseccion.biseccion(a, b, TOL, 100)
  2. Biseccion.dibujar(a,b,TOL, 100)

Bien, en estas dos llamadas lo que hacemos es llamar a las funciones bisección y dibujar residentes en nuestro fichero Bisección.py, en ambas llamadas pasamos 3 parámetros variables: a, b y TOL; y uno fijo, el 100, que va a ser el número máximo de iteraciones que le vamos a permitir que realice la Bisección sin que entre en un bucle quasi-infinito, si queréis podéis modificar éste número aumentándolo o disminuyendo su cantidad según os apetezca. Como el método de la bisección es muy lento, si no el más lento de todos, lo mejor será dejarlo entre 100 y 150, para Newton se puede bajar mucho más. la rapidez de la bisección va a depender muchísimo de la longitud del intervalo, cuanto mayor sea esta longitud mayor será el número de iteraciones a realizar, ya que lo que hacemos es dividir los diferentes subintervalos por la mitad.

Bien, pasemos al código del fichero evaluar.py, que es el siguiente:

  1. from math import *
  2.  
  3. lista_segura = ['math','acos', 'asin', 'atan', 'atan2', 'ceil', 'cos', 'cosh', 'de grees', 'e', 'exp', 'fabs', 'floor', 'fmod', 'frexp', 'hypot', 'ldexp', 'log', 'log10', 'modf', 'pi', 'pow', 'radians', 'sin', 'sinh', 'sqrt', 'tan', 'tanh']
  4. dicc_seguro = dict([ (k, locals().get(k, None)) for k in lista_segura ])
  5. funcion = raw_input("Introduce la funcion\n")

Como véis es muy corto, consta de una importación de la librería math de Python, de la declaración de una lista formada por los nombres de las principales funciones matemáticas utilizadas en cualquier lenguaje de programación, de un diccionario formado a partir de la lista y con el cual vamos a poder hacer también ciertos ajustes de variable. Finalmente esta la línea de toma de datos de la función matemática a considerar. Todo ello nos va a venir de perlas para poder utilizar eval en el fichero Bisección.py de nuestro proyecto, el cual aquí os presento ya su código:

  1. #!/usr/bin/python
  2. import evaluar
  3. from pylab import *
  4. from Numeric import *
  5.  
  6. def biseccion(a, b, TOL, N):
  7.     evaluar.dicc_seguro['x']=a
  8.     fa = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro)
  9.  
  10.     vectorx = zeros(N, Float64)
  11.     vectory = zeros(N, Float64)
  12.  
  13.     i = 1
  14.     while i<=N :
  15.         p = (a+b)/2.0
  16.         vectorx[i-1] = p
  17.         evaluar.dicc_seguro['x']=p
  18.         fp = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro)
  19.         vectory[i-1]=fp
  20.         if (fp == 0.0) or ((b-a)/2.0)<tol:
  21.             print "La raiz buscada es: ",p, "con", i-1, "iteraciones"
  22.             break
  23.         i = i+1
  24.         if (fa*fp)>0 :
  25.             a = p
  26.         else :
  27.             b = p
  28.  
  29.     return [vectorx, vectory]
  30.  
  31. def dibujar(a,b, TOL, N):
  32.   x = arange(a,b,0.1)
  33.   vectores=biseccion(a, b, TOL, N)
  34.  
  35.   subplot(211)
  36.   plot(x, eval(evaluar.funcion), linewidth=1.0)
  37.   xlabel('Abcisa')
  38.   ylabel('Ordenada')
  39.   title('Metodo Biseccion con f(x)=' + evaluar.funcion)
  40.   grid(True)
  41.   axhline(linewidth=1, color='r')
  42.   axvline(linewidth=1, color='r')
  43.  
  44.   subplot(212)
  45.   plot(vectores[0], vectores[1], 'k.')
  46.   xlabel('Abcisa')
  47.   ylabel('Ordenada')
  48.   grid(True)
  49.   axhline(linewidth=1, color='r')
  50.   axvline(linewidth=1, color='r')
  51.  
  52.   show()

Lo primero a comentar es que en las primeras líneas importamos nuestro fichero evaluar.py la librería pylab para poder dibujar nuestras funciones. Después vemos que hemos declarado dos funciones, una llamada bisección y otra llamada gráfica, para ello en Python anteponemos la palabra clave def, seguida del nombre queremos ponerle a la función y de la lista de parámetros cerrados entre paréntesis, en ambas funciones que hemos declarado pasamos cuatro parámetros: a, b, TOL y N.
De la función bisección no voy a entrar en detalles de cómo funciona, sólo decir que lo que hace es dividir el intervalo de partida en subintervalos partidos por el punto medio, quedándose con aquel subintervalo que cumpla el Teorema de Bolzano. Obviamente, antes de hacer nada hemos de escoger funciones que sean continuas en el intervalo de partida, ya que es una de las hipótesis de Bolzano, la otra es la del cambio de signo entre las imágenes de los extremos de los diferentes subintervalos.
Veamos qué significan las siguientes dos líneas en términos de Python:

  1. evaluar.dicc_seguro['x']=a
  2. fa = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro)

Bien, en la primera línea llamamos a nuestro diccionario del fichero
evaluar.py, el diccionario le hemos llamado dicc_seguro, y lo que hacemos es que la variable ‘x’, que se asigna por defecto a nuestra función, pase a ser identificada por ‘a’, para así poder calcular f(a) en la siguiente línea mediante la función eval, que tiene una estructura extraña, el primer parámetro, evaluar.funcion, lo que hace es tomar el valor de funcion de evaluar.py (por ejemplo podría ser cos(x)-x ), el segundo lo copiáis como esta porque no lo he entendido muy bien, pero parece ser que es por motivos de seguridad, y el tercero lo que hacemos es llamar a nuestro diccionario con el valor de nuestra variable ‘a’. Esto lo hacemos también para calcular ‘fp’.
Otras dos líneas a explicar son:

  1. vectorx = zeros(N, Float64)
  2. vectory = zeros(N, Float64)

Aquí declaramos dos vectores inicializados sus elementos con valor 0.0, de carácter Float64, uno va a ser las coordenadas x y el otro las coordenadas y, que especificamos más adelante como:

  1. vectorx[i-1] = p
  2. vector[i-1] = fp

Bien, vectorx y vectory van a ser los puntos aproximativos del método de la bisección y que vamos a utilizar lo que he llamado la nube de puntos aproximativa de la bisección mediante la función grafica. Estos vectores no son necesarios para la bisección, pero sí me hacen falta para la nube de puntos. De la función bisección falta comentar la siguiente línea:

  1. return [vectorx, vectory]

Esto es lo que más me ha costado de hacer en el código, y que gracias a mi amigo Kawalero he podido resolver, bueno me lo ha resuelto él. Y es que no sabía cómo devolver más de un dato de una función, y la solución es mediante una lista, en este caso la lista esta formada por dos elementos que son dos vectores. y que los devolvemos con la instrucción return, muy conocida en C/C++. Ahora está el cómo llamar a vectorx y vectory en la función grafica, lo cual resolvemos así:

  1. vectores=biseccion(a, b, TOL, N)
  2. plot(vectores[0], vectores[1], 'k.')

Montamos una lista llamada ‘vectores’, mediante la llamada a la función bisección, y luego vectorx será vectores[0] y vectory será vectores[1], los cuales utilizamos en la función plot para dibujar la nube de puntos aproximativa.
Sobre la función grafica no voy a explicar mucho más, si estáis acostumbrados/as al uso de Matlab, Octave o Scilab no os resultará difícil su comprensión, aun así os dejo el enlace Pylab

Y bueno, podéis ver el código en ejecución con este ejemplo:

Introduce la funcion
sin(x)
Introduce el extremo inferior del intervalo
-0.8
Introduce el extremo superior del intervalo
0.2
Introduce la tolerancia del metodo
0.0001
La raiz buscada es: -1.220703125e-05 con 13 iteraciones

Pantallazo-Figure 1 

Como podéis ver obtenemos un acabado muy completo de la función, lo he dividido en dos partes, además por defecto podemos guardar la imagen en el formato que queramos, así como escalarla, hacer zoom y etc.

Cosas Que Me Quedan Por Hacer

  • Implementar en python el código de los demás métodos de aproximación.
  • Sustituir la recogida de datos por terminal por código basado en wxWidgets.
  • Bien, podéis consultar si queréis ver la bisección en forma recursiva en ForosDelWeb, en donde pregunté algunos lios que me hice, también podréis ver otra forma de devolver más de un dato con clases, gracias a Kawalero. A todos los que me ayudaron muchas gracias

    Ala saludotes, me voy a ver si me refresco que hace mucho calor peña :-h

    9 pensamientos en “La Bisección Con Python

    1. Sanatas

      Pues si que parece que te gusta entrar rápido en faena, peazo de post te has currado, enhorabuena.

    2. Cristobal

      Jejejeje, como dicen en mi terreta “pensado y hecho”. Muchas gracias tio A ver si no aprieta tanto la calor y voy escribiendo más cosas de este tipo.

    3. ANdres

      Hola buenas!
      Primero buena currada!!
      Queria preguntar una cosilla: El modulo Biseccion, ¿dónde lo tienes guardado? ¿Esta donde el resto de las librerías?
      Muchas gracias y Un saludo!!

    4. Cristobal

      Hola Andrés, gachias por el cumplido
      Lo del módulo de la bisección lo tengo todo guardado en el mismo directorio/carpeta, el proyecto lo hice con Eclipse utilizando la extensión para Python.
      De nadas

    5. ANdres

      Hola!!
      Yo es que tengo también varios modulos en un mismo directorio y no consigo importar unos en otros. Trabajo con Windows XP…
      Muchas gracias (y gracias por tu anterior respuesta)

    6. ANdres

      Hola de nuevo!!
      Parece que ya me deja importar. No debía tener el programa principal en el mismo directorio que los modulos que quería usar…
      Muchas gracias, seguiremos currelando.

      Saludos!!!

    7. pedaso_de_loco

      epale Cristobal… ya que me da la impresion de que sabes mucho Python.. yo ahora es que voy a trabajar en ese lenguaje de programacion, me gustaria.. si esta a su alcanc que me eche una mano… y aprender!

    8. Pingback: Métodos De Aproximación De Raíces Con Python

    9. german

      no puedo ejecutarlo, por q me dice q no poseo el modulo numeric, pero si lo tengo instalado, y lo pruebo aparte y funciona le modulo. por q sera?

    Los comentarios están cerrados.