Archivos de Categoría: Python - Paginas 2

Trazadores Cúbicos De Tipo 1 Con Python

Hola, en general un trazador o función spline en un intervalo [a,b], en matemáticas, son polinomios que están definidos en subintervalos de [a,b] que cumplen ciertas condiciones de regularidad; normalmente cumplen lo que se llama ser función de clase C^n([a,b]) , es decir, que son continuas tanto la función como sus derivadas hasta orden n en [a,b], lo cual es bastante obvio en polinomios.

Los trazadores se utilizan cuando por ejemplo tenemos un conjunto de puntos y queremos encontrar la función que pasa por todos esos puntos y que mejor se aproxima, dicha función a calcular es lo que llamamos trazador o función spline. Esta técnica de los trazadores es muy útil en informática para dibujar figuras no regulares. Por ejemplo, se puede utilizar para dibujar un perro o un edificio con tan sólo dar unos cuantos puntos representativos de dichas figuras y buscando su trazador.

El uso de trazadores en informática es muy útil ya que nos permite encontrar funciones que pasan por una cantidad finita de puntos mediante funciones a trozos, todas ellas polinómicas de grado pequeño, con lo cual el cómputo de operaciones de los algoritmos se reduce drásticamente, además de que dichas funciones a trozos (los trazadores) son funciones “suaves”.

Otro campo matemático parecido a los splines son los NURBS, los cuales nos permiten representar figuras complicadas y libres mediante interacción con el usuario mediante el movimiento de los nodos. Son muy utilizados en campos como la animación, y también se les conoce como B-splines, por ser un caso particular de los que veremos aquí. Si habéis programado con OpenGL sabréis de lo que estoy hablando.
NURB

Pero dejemos de lado estas profundas disquisiciones y profundicemos un poco en los trazadores cúbicos, los cuales son un buen comienzo para adentrarse en este apasionante mundo de la aproximación de figuras libres utilizando las matemáticas. Empecemos definiendo que es una función spline en general.

Sea \Delta=\{ a=x_0 <  x_1   < .... <  x_n= b \} una partición del intervalo [a,b]. S_{\Delta} : [a,b]\longrightarrow\mathbb{R} es una función spline de orden k\in\mathbb{N} asociada a \Delta si S_{\Delta}\in C^{k-1} ([a,b]) y si S_{\Delta} coincide en cada intervalo [ x_i , x_{i+1} ] , i=0,1,...,n-1, con un polinomio de grado \le  k .
Notemos que para k=1 se obtiene una función lineal a trozos, que todos llamamos poligonal.
En éste artículo me centraré sólo en el caso k=3, conocidos como trazadores cúbicos o splines cúbicas. Así pues nuestras splines serán de clase C^2 y de grado menor o igual que 3. Y en concreto a las de tipo I que definiré en breve.
Dado y = ( y_0 , y_1 , …. y_n ) \in\mathbb{R}^{n+1} denotaremos por S_\Delta (y, .) a una función spline cúbica de interpolación que verifica S_\Delta (y,x_i ) = y_i para i = 0,1,…,n-1.
Hay que recalcar que estas condiciones no determinan de forma única la spline cúbica. Si os estáis preguntando qué pinta el y es muy sencillo, en nuestros problemas partiremos de una serie de puntos concretos, por ejemplo, consideremos que queremos interpolar los puntos (0,1), (1,2), (4,0); pues bien, cada abscisa de los puntos representan la partición y las ordenadas de cada punto representan la y, que lógicamente estos son los extremos de cada subintervalo de cada función a trozo de que se compone nuestra función spline.
Para determinar únivocamente las funciones spline se imponen, en general, 3 tipos de condiciones. Como en éste artículo sólo me voy a centrar en las de tipo I, su condición de unicidad va a ser:
S^{\prime\prime}_\Delta (y,a) = S^{\prime\prime}_\Delta (y,b) = 0. Se suelen denominar condiciones naturales.
Ya hablaré de los otros dos tipos en otra ocasión, a modo de información os diré que las de tipo III se utilizan para interpolar funciones periódicas de periodo (b-a).
Para calcular explícitamente los trazadores definiremos primero los momentos de S_{\Delta}(y,.) mediante un vector M cuyas componentes vienen expresadas como M_j = S^{\prime\prime}_{\Delta}(y,x_{j}) para j=0,1,…,n-1.
Denotemos h_j = x_{j+1}-x_j j=0,1,…,n-1 Para calcular los momentos utilizaremos loa siguiente expresión:
\mu_j M_{j-1} + 2M_j + \lambda_j M_{j+1} = d_j  1\le j \le n-1
Donde los vectores \mu,\lambda,d vienen definidos por las expresiones:
\lambda_j = \frac{h_{j+1}}{h_j + h_{j+1}}
\mu_j = 1-\lambda_j
d_j = \frac{6}{h_j + h_{j+1}} \cdot ( \frac{y_{j+1}-y_j}{h_{j+1}}- \frac{y_j-y_{j-1} }{h_j} )
para 1 \le j \le n-1
Para los extremos (j=0, j=n) imponemos las condiciones de tipo I S^{\prime\prime}_\Delta (y,a) =  S^{\prime\prime}_\Delta(y,b) = 0 y obtenemos :
M_0 = M_n = 0
\lambda_0 = d_0 = \mu_n = d_n = 0
Si reunimos todas las fórmulas obtenidas tenemos que llegamos a un sistema lineal AM = d donde A es una matriz tridiagonal en la que su diagonal principal esta formada todos su elementos por el 2, su diagonal superior lo forman el vector \lambda y su diagonal inferior esta formada por el vector \mu .
Con todos estos datos ya podemos pasar a realizar un programa en Python sobre un caso concreto para ver cómo funciona. Las demostraciones de cómo se obtienen las fórmulas anteriores las podéis obtener en cualquier libro de texto universitario que se dedique a la rama matemática de Aproximación Numérica bajo el apartado de interpolación polinómica bajo splines o trazadores.
EJEMPLO
Vamos a determinar la función spline cúbica de tipo I que interpola los puntos:
(0,1), (1,2),(3,0)
Para ello he adaptado el método de trazadores cúbicos de tipo I a un sencillo programa escrito bajo lenguaje Python utilizando las librerías gráficas Matplotlib para poder dibujar la función spline. Su código es el siguiente:

from Numeric import *
from math import *
from numpy.linalg import *
from pylab import *

vectorx = array((0.,1.,3.), Float)
vectory = array((1.,2.,0.), Float)
h = array((0.,0.), Float)
lamda = array((0.0, 0.0),Float)
mu = array((0.0, 0.0),Float)
d = array((0.0, 0.0,0.0),Float)
A = array(([0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]), Float)

i = 0
n = 0

while i<len(h) :
	h[i] = vectorx[i+1]-vectorx[i]
	i=i+1
i=1
while i<len(h) :
	lamda[i]=h[i]/(h[i-1]+h[i])
	i=i+1

i = 1
while i<len(h) :
	mu[i-1]=1.-lamda[i]
	i=i+1

i=1
while i<len(h) :
	d[i] = (6./(h[i-1]+h[i]))*((vectory[i+1]-vectory[i])/h[i]-(vectory[i]-vectory[i-1])/h[i-1])
	i=i+1

i=0
n=1
while i<len(mu) :
	A[i][n]=lamda[n-1]
	i=i+1
	n=n+1

i=1
n=0
while n<len(mu) :
	A[i][n]=mu[i-1]
	i=i+1
	n=n+1

i=0
while i<=len(mu) :
	A[i][i]=2.
	i=i+1

M = solve(A,d)

def funcionS1() :
	x = arange(0,3,0.1)
	m = len(x)
	s1 = array(([0.,0.]),Float)
	s2 = array(([0.0,0.]),Float)
	S = zeros(m, Float64)
	i = 0
	l = 0
	while i<=1 :
		while l<m :
			s1[i]=(vectory[i+1]-vectory[i])/(h[i+1])-(1./6.)*(2*M[i]+M[i+1])*h[i+1]
			s2[i]=(1./6.)*((M[i+1]-M[i])/h[i+1])
			S[l] = vectory[i]+s1[i]*(x[l]-vectorx[i])+0.5*M[i]*pow(x[l]-vectorx[i],2.)+s2[i]*pow(x[l]-vectorx[i],3.)
			l = l+1
		i = i+1
        plot(x, S, linewidth=1.0)
        xlabel('Abcisa')
        ylabel('Ordenada')
        title('Metodo I De Los Splines ')
        grid(True)
        axhline(linewidth=1, color='r')
        axvline(linewidth=1, color='r')
        show()

funcionS1()

El código del programa no tiene mayor misterio, a destacar la instrucción para resolver el sistema lineal AM=d, que es:

M = solve(A,d)

He utilizado la función solve de las librerías Numpy, en concreto, las específicas al álgebra lineal, y que para importarlas a nuestros programas se debe hacer de una forma especial al inicio de nuestro código mediante la instrucción:

from numpy.linalg import *

He decidido utilizar la función solve porque nuestro sistema lineal es muy pequeño, pero en casos grandes es aconsejable utilizar métodos de cálculo de sistemas lineales con matrices de coeficientes tridiagonales, para así disminuir drásticamente el número de operaciones que realice nuestro código.
De igual forma cuando calculamos los puntos interpoladores de nuestra función spline, si la cantidad de puntos a calcular es grande, es aconsejable utilizar el Algoritmo de Horner para su cálculo, lo cual es posible gracias a que nuestras funciones spline son funciones polinómicas definidas a trozos.
El resultado del programa lo podéis ver en la siguiente imagen:
Trazador Cúbico Tipo I

Espero pronto poder escribir sobre los trazadores cúbicos de tipo II y III

Saludos :-h

Métodos De Aproximación De Raíces Con Python

Hola, siguiendo con el artículo de La Bisección con Python he añadido dos nuevos métodos para aproximar raíces de funciones. He añadido los métodos del Punto Fijo y el de la Secante. Como el código lo que tiene de nuevo respecto al anterior es aprender en Python a hacer un menú me he decidido a comentarlo un poco.
Quisiera dejar claro que éste artículo no va a tratar sobre una explicación de estos métodos iterativos con profundidad, ya que considero que hay mucho material escrito sobre ello y se puede encontrar fácilmente en cualquier biblioteca pública universitaria de ciencias, en librerías y en artículos gratuitos en formato pdf en la red. Además dichos métodos no contienen ningún misterio.
Bien, pasemos a hablar un poco del código que he realizado, el cual se compone de 6 ficheros con extensión .py. El primero es el archivo principal del proyecto, el cual he llamado Metodos.py y cuyo código es:

#!/usr/bin/python
# coding: latin-1

#import os, sys

import Menu
from Numeric import *

Menu.elmenu()

Este código no tiene ningún misterio, lo único que hace es cargar python con codificación de latin. Nos importa el fichero Menu.py que hemos implementado nosotros mismos, y de Menu.py llamamos a la función elmenu() con la instrucción Menu.elmenu()
Pasemos a ver el código del fichero Menu.py:

#!/usr/bin/python
# coding: latin-1

from math import *

import Biseccion
import PuntoFijo
import Secante

def op1 ():
    a = float(input("Introduce el extremo inferior del intervalo\n"))
    b = float(input("Introduce el extremo superior del intervalo\n"))
    TOL = float(input("Introduce la tolerancia del metodo\n"))
    Biseccion.biseccion(a, b, TOL, 100)
    Biseccion.dibujar(a, b, TOL, 100)
def op2 ():
    po = float(input("Introduce el valor inicial\n"))
    TOL = float(input("Introduce la tolerancia del metodo\n"))
    PuntoFijo.puntofijo(po, TOL, 100)
    PuntoFijo.dibujar(po, TOL, 100)
def op3 ():
    po = float(input("Introduce el primer valor inicial\n"))
    p1 = float(input("Introduce el segundo valor inicial\n"))
    TOL = float(input("Introduce la tolerancia del metodo\n"))
    Secante.secante(po, p1, TOL, 100)
    Secante.dibujar(po, p1, TOL, 100)
def errhandler ():
   print "Tu eleccion no ha sido la correcta\n"

#menu es un diccionario que nos ayudara para hacer un menu por terminal
def elmenu ():
    menu = {"1": op1, "2": op2, "3": op3}
    opcion = raw_input("Por favor introduce la opcion numerica que desees\n")
    menu.get(opcion, errhandler)()

Lo primero que hacemos es importar la librería matemática “math” de Python y los métodos de la Bisección, Punto Fijo y Secante, los cuales los he escrito en tres ficheros .py diferentes, lo hacemos con el trozo de código:

from math import *

import Biseccion
import PuntoFijo
import Secante

Seguidamente pasamos a definir las cosas que realizarán cada una de las opciones de que se compone nuestro menú, por ejemplo, con la opción 1 pediremos por terminal los valores de a, b y TOL, además mandamos que se ejecuten las funciones bisección y dibujar que se encuentran en Biseccion.py, el trozo de código correspondiente es:

def op1 ():
    a = float(input("Introduce el extremo inferior del intervalo\n"))
    b = float(input("Introduce el extremo superior del intervalo\n"))
    TOL = float(input("Introduce la tolerancia del metodo\n"))
    Biseccion.biseccion(a, b, TOL, 100)
    Biseccion.dibujar(a, b, TOL, 100)

De igual forma se procede con op2() y op3(), sólo hay que tener en cuenta de anteponer la palabra clave def. He añadido una opción de salida del menú en caso de error:

def errhandler ():
   print "Tu eleccion no ha sido la correcta\n"

Finalmente definimos la función elmenu(), la cual es la que hemos llamado en Metodos.py. Esta función se compone de un diccionario llamado menu, así pues, un menú en Python lo podemos definir mediante un diccionario como el que se ha definido en este caso. Luego la función se compone de una toma de las distintas opciones del menú mediante raw_input y menu.get. El trozo de código es:

def elmenu ():
    menu = {"1": op1, "2": op2, "3": op3}
    opcion = raw_input("Por favor introduce la opcion numerica que desees\n")
    menu.get(opcion, errhandler)()

Seguidamente viene el código de Biseccion.py, el cual es:

#!/usr/bin/python

import Evaluar
from pylab import *
from Numeric import *

def biseccion(a, b, TOL, N):
    Evaluar.dicc_seguro['x']=a
    fa = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
    vectorx = zeros(N, Float64)
    vectory = zeros(N, Float64)

    i = 1
    while i<=N :
        p = (a+b)/2.0
        vectorx[i-1] = p
        Evaluar.dicc_seguro['x']=p
        fp = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
        vectory[i-1]=fp
        if (fp == 0.0) or ( (b-a)/2.0)<tol:
            break
        i = i+1
        if (fa*fp)>0 :
            a = p
        else :
            b = p
    print "La raiz buscada es: ",p, "con", i-1,"iteraciones"
    return [vectorx, vectory]

def dibujar(a,b,TOL, N):
  x = arange(a,b,0.1)
  vectores=biseccion(a, b, TOL, N)

  subplot(211)
  plot(x, eval(Evaluar.funcion), linewidth=1.0)
  xlabel('Abcisa')
  ylabel('Ordenada')
  title('Metodo Biseccion con f(x)=' + Evaluar.funcion)
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  subplot(212)
  plot(vectores[0], vectores[1], 'k.')
  xlabel('Abcisa')
  ylabel('Ordenada')
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  show()

Este código ya lo expliqué en el anterior artículo, se compone de dos funciones, en una de ellas realizamos el método de la bisección y en la otra función dibujamos la gráfica de la función y sus aproximaciones.
De igual forma he escrito el método del punto fijo con PuntoFijo.py

#!/usr/bin/python

import Evaluar
from pylab import *
from Numeric import *

def puntofijo(po,TOL, N):
    vectorx = zeros(N, Float64)
    vectory = zeros(N, Float64)

    i = 1
    while i<=N :
        vectorx[i-1] = po
        Evaluar.dicc_seguro['x']=po
        fp = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
        vectory[i-1]=fp
        if fabs(po-fp)<tol:
            print "La raiz buscada es: ",po, "con", i-1, "iteraciones"
            break
        i = i+1
        po = fp
    return [vectorx, vectory]

def dibujar(po,TOL, N):
  x = arange(po-2,po+2,0.1)
  vectores=puntofijo(po, TOL, N)

  subplot(211)
  plot(x, eval(Evaluar.funcion), linewidth=1.0)
  xlabel('Abcisa')
  ylabel('Ordenada')
  title('Metodo Punto Fijo con f(x)= x - ' + Evaluar.funcion)
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  subplot(212)
  plot(vectores[0], vectores[1], 'k.')
  xlabel('Abcisa')
  ylabel('Ordenada')
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  show()

Este código hace lo mismo que Biseccion.py pero con el método del punto fijo. Para el de la secante lo he puesto en Secante.py y su código es:

#!/usr/bin/python

import Evaluar
from pylab import *
from Numeric import *

def secante(po, p1, TOL, N):
    Evaluar.dicc_seguro['x']=po
    qo = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
    Evaluar.dicc_seguro['x']=p1
    q1 = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
    vectorx = zeros(N, Float64)
    vectorx1 = zeros(N, Float64)
    vectory = zeros(N, Float64)
    vectory1 = zeros(N, Float64)

    i = 2
    while i<=N :
        p = p1-(q1*(p1-po)/(q1-qo))
        vectorx[i-2] = po
        vectorx1[i-2] = p1
        Evaluar.dicc_seguro['x']=p
        fp = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
        vectory[i-2]=qo
        vectory1[i-2]=q1
        if fabs(po-p1)<tol :
        	print "La raiz buscada es: ",p, "con", i-2,"iteraciones"
        	break

       	i = i+1
       	po = p1
       	qo = q1
       	p1 = p
       	q1 = fp

    return [vectorx, vectory, vectorx1, vectory1]

def dibujar(po, p1, TOL, N):
  x = arange(po-2,po+2,0.1)
  vectores=secante(po, p1, TOL, N)

  subplot(211)
  plot(x, eval(Evaluar.funcion), linewidth=1.0)
  xlabel('Abcisa')
  ylabel('Ordenada')
  title('Metodo Secante con f(x)= ' + Evaluar.funcion)
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  subplot(212)
  plot(vectores[0], vectores[1], 'r.', vectores[2], vectores[3], 'b.')
  xlabel('Abcisa')
  ylabel('Ordenada')
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  show()

Y lo que hace el código pues idem de idem.
Finalmente queda el código de Evaluar.py, que es:

from math import *

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']
dicc_seguro = dict([ (k, locals().get(k, None)) for k in lista_segura ])

print"METODOS ITERATIVOS PARA APROXIMAR RAICES."
print "________________________________________\n"
print "1. Metodo de la Biseccion\n"
print "2. Metodo del Punto Fijo\n"
print "3. Metodo de la Secante\n"

funcion = raw_input("Introduce la funcion para la cual vas a aplicar uno de los metodos\n")

Con dicho código establec emos primero una lista de las principales funciones y constantes matemáticas para que el usuario cuando las introduzca por terminal puedan ser interpretadas correctamente por nuestro programa, es decir, que si nuestro usuario introduce por ejemplo:
\frac{cos(x)}{x}
Pues que lo interprete nuestro código como la fracción entre cos(x) y x, cuando le pasamos ello a la variable funcion con raqw_input.
Bueno os dejo un pantallazo de cómo quedaría el resultado del programa utilizando la función cos(x)-x con el método de la secante, usando como puntos iniciales p_0 = 0.5, p_1 = \pi/4

Programa Python 

Saludos :-h

Ampliando Mi Documentación De Python

Hola, he ampliado mi documentación de Python (¡ gracias Sanatas! ) y he decidido compartirla con vosotros, el otro día ya puse un enlace de descarga, os pongo ahora el nuevo enlace de descarga.

DESCARGAR DOCUMENTACIÓN DE PYTHON

A parte de apuntes y tutoriales he añadido ejemplos hechos con librerías en Python que podemos instalar desde Synaptic pero que luego es engorroso acceder a ellos, por ejemplo, los ejemplos de wx2.8Python hay además que ejecutar un script de descompresión en modo superusuario para tener acceso a ellos. Os detallo un poco lo que contiene la documentación.

APUNTES TEÓRICOS

  1. Python para todos de Zootropo
  2. Introducción a la programación con Python de la Univ. Jaume I de Castellón
  3. Documentación oficial de Python 2.5.2 (lo tenéis en una carpeta)
  4. Lo anterior en español y en formato HTML
  5. Apuntes de PyGTK+2 en inglés y castellano
  6. Tutorial de PyGTK+2 en formato PDF.
  7. Documentación científica de Python en PDF
  8. Documentación de Python con KDE4
  9. Tutorial de PyQT de ELCODIGOK
  10. Tutorial de las librerías NumPY en formato PDF
  11. Tutorial Dive Into Python en PDF, muy extenso, y con ejercicios propuestos.
  12. Tutorial de las librerías científicas SciPY
  13. Tutorial de las librerías gráfico-matemáticas Matplotlib
  14. Tutorial de las librerías gráficas Python-Visual
  15. Tutorial de las librerías TKinter en PDF
  16. Libro Python In Action sobre wxPython en PDF
  17. Y una peequeña licencia que me he permitido de incluir porque lo tengo suelto sobre Métodos Numéricos con FORTRAN

EJEMPLOS PYTHON

  1. De wxPython en la versión 2.8
  2. De las librerías Matplotlib
  3. De las librerías gráficas Mayavi
  4. Los programas de Dive Into Python
  5. Colección de ejemplos de iniciación en Python, tenéis un montón (bajo la carpeta programas-python)
  6. Ejemplos de PyGTK2+ (esán dentro de la carpeta pygtk2tutorial-es/examples)
  7. Ejemplos de Python-Visual

Y bueno esta es mi colección de apuntes sobre Python. Espero que os sean de ayuda para aprender peña ;)

Saludos :-h

Mis Apuntes De Python

Hola, como ya muchos sabéis este verano me ha dado por aprender Python, así que me he recopilado unos cuantos apuntes en formato PDF sobre Python y me los he subido a mi cuenta de MediaFire para tenerlos a mano por si los pierdo. He pensado que a alguien le puede venir bien así que los podéis descargar del siguiente enlace:

DESCARGAR APUNTES PYTHON

El fichero comprimido contiene 9 ficheros, entre ellos un libro de cómo aprender ha hacer widgets con Python en inglés pero muy bueno, también está el de Zootropo que es un manual muy bueno sobre iniciación en Python. Hay tres tutoriales sobre las librerías Numpy, SciPy y MatPlotLib. También hay un tutorial sobre diferentes consideraciones para realizar proyectos Python en ciencias, luego está la guía de Python2.5 y otro que es un libro sobre Python. Si a alguien le sirve pues guay ;)

Saludos :-h

La Bisección Con Python

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:

#!/usr/bin/python
# coding: latin-1

import os, sys
import Biseccion
from Numeric import *

a = float(input("Introduce el extremo inferior del intervalo\n"))
b = float(input("Introduce el extremo superior del intervalo\n"))
TOL = float(input("Introduce la tolerancia del metodo\n"))
Biseccion.biseccion(a, b, TOL, 100)
Biseccion.dibujar(a,b,TOL, 100)

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

#!/usr/bin/python
# coding: latin-1
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:

import Biseccion
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:

a = float(input("Introduce el extremo inferior del intervalo\n"))
b = float(input("Introduce el extremo superior del intervalo\n"))
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:

Biseccion.biseccion(a, b, TOL, 100)
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:

from math import *

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']
dicc_seguro = dict([ (k, locals().get(k, None)) for k in lista_segura ])
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:

#!/usr/bin/python
import evaluar
from pylab import *
from Numeric import *

def biseccion(a, b, TOL, N):
    evaluar.dicc_seguro['x']=a
    fa = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro)

    vectorx = zeros(N, Float64)
    vectory = zeros(N, Float64)

    i = 1
    while i<=N :
        p = (a+b)/2.0
        vectorx[i-1] = p
        evaluar.dicc_seguro['x']=p
        fp = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro)
        vectory[i-1]=fp
        if (fp == 0.0) or ((b-a)/2.0)<tol:
            print "La raiz buscada es: ",p, "con", i-1, "iteraciones"
            break
        i = i+1
        if (fa*fp)>0 :
            a = p
        else :
            b = p

    return [vectorx, vectory]

def dibujar(a,b, TOL, N):
  x = arange(a,b,0.1)
  vectores=biseccion(a, b, TOL, N)

  subplot(211)
  plot(x, eval(evaluar.funcion), linewidth=1.0)
  xlabel('Abcisa')
  ylabel('Ordenada')
  title('Metodo Biseccion con f(x)=' + evaluar.funcion)
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  subplot(212)
  plot(vectores[0], vectores[1], 'k.')
  xlabel('Abcisa')
  ylabel('Ordenada')
  grid(True)
  axhline(linewidth=1, color='r')
  axvline(linewidth=1, color='r')

  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:

evaluar.dicc_seguro['x']=a
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:

vectorx = zeros(N, Float64)
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:

vectorx[i-1] = p
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:

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í:

vectores=biseccion(a, b, TOL, N)
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

    Page 2 of 212