Duda con array

[tofol]

Hola estoy haciendo un programa que me calcule raices con el método de la bisección, y quiero que me vaya guardando los valores de p y fp en un vector utilizando array de Numpy. El programa no me echa ningún error, pero en p no me guarda los valores que debería guardar.
he dividido el programa en tres ficheros:

Fichero evaluar.py

Código:

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")

Fichero Metodos.py (este es el fichero main del programa)

Código:

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

import os, sys
import Biseccion

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"))
N = int(input("Introduce el numero maximo de iteraciones a realizar\n"))
Biseccion.biseccion(a, b, TOL, N)
Biseccion.dibujar(a,b)

El fichero Biseccion.py

Código:

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

def biseccion(a, b, TOL, N):
    evaluar.dicc_seguro['x']=a
    fa = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro) 
    vectorx = zeros(N)
    i = 1
    while i<=N :
        p = (a+b)/2
        #vectorx[i] = p
        print p
        for i in range(1, N):
            vectorx[i]=p    
        evaluar.dicc_seguro['x']=p
        fp = eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro)
        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
    print "Vector x",vectorx
    
def dibujar(a,b):
  t = arange(a,b,0.1)
  evaluar.dicc_seguro['x']=t
  plot(t, eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro), linewidth=1.0)
  xlabel('Abcisa')
  ylabel('Ordenada')
  title('Metodo Biseccion con f(x)=' + evaluar.funcion)
  grid(True)
  axhline(linewidth=2, color='r')
  axvline(linewidth=2, color='r')
  show()

El problema reside en este último fichero en la función biseccion, he optado por poner el bucle este

Código:

 for i in range(1, N):
            vectorx[i]=p

Pero no me guarda los diferentes valores de p, sólo me guarda un único valor para todos los valores de vectorx

Pongo aquí un ejemplo de salida del programa:

Introduce la funcion
x-1.3333
Introduce el extremo inferior del intervalo
1
Introduce el extremo superior del intervalo
2
Introduce la tolerancia del metodo
0.00001
Introduce el numero maximo de iteraciones a realizar
10
1.5
1.25
1.375
1.3125
1.34375
1.328125
1.3359375
1.33203125
1.333984375
1.3330078125
1.33349609375
1.33325195312
1.33337402344
1.33331298828
1.3332824707
1.33329772949
1.33330535889
La raiz buscada es: 1.33330535889
Vector x [ 0. 1.33330536 1.33330536 1.33330536 1.33330536 1.33330536
1.33330536 1.33330536 1.33330536 1.33330536]

Como se puede apreciar en Vector x debería de tener los valores 1.5, 1.25, 1.375 etcétera, pero no es así. ¿Alguien me puede decir qué estoy haciendo mal?
Aclarar que el método de la bisección me va correcto pero necesito vectorx para poder dibujar los puntos de aproximación del método en una gráfica.
Por cierto, esto no son mis deberes de clase, estoy aprendiendo Python por mi cuenta y ahora no se a quién preguntar porque mis tiempos de facultad ya hace años que acabaron y no eran de programador

Muchas gracias

[vermells]

hola tofol yo tambien uso python auqnue no se mucho mucho,lo que veo en tu codigo es que declaras vectorx= zero(N) pero no tendria que ser vectorx=[] primero?
porque sino no hay array,bueno no se si es eso,un saludo chao

[vermells]

bueno me precipite y ese no es el problema,la verdad es que no se cual es,haber si alguien da la solucion un saludo

[tofol]

Hola vermells, gracias por tu ayuda, ya lo he resuelto, era un problema de índices. Pongo aquí el Biseccion.py corregido

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 = vectorx
    i = 1
    while i<=N :
        p = (a+b)/2
        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
    print "Vector x = ", vectorx
    print "Vector y = ", vectory
    
def dibujar(a,b):
  t = arange(a,b,0.1)
  evaluar.dicc_seguro['x']=t
  plot(t, eval(evaluar.funcion, {"__builtins__":None}, evaluar.dicc_seguro), linewidth=1.0)
  xlabel('Abcisa')
  ylabel('Ordenada')
  title('Metodo Biseccion con f(x)=' + evaluar.funcion)
  grid(True)
  axhline(linewidth=2, color='r')
  axvline(linewidth=2, color='r')
  show()

Sólo me falta una cosa y es poder utilizar vectorx y vectory dentro de la función dibujar de Biseccion.py, como vectorx y vectory sólo los puedo calcular en la función biseccion, no se cómo llamarlos desde dibujar. Lo he intentado con return multiply pero nada de nada. ¿Alguien sabe cómo hacerlo? Por si acaso os explico, quiero pillar vectorx, vectory en dibujar para ponerlos en la instrucción plot, y así conseguir que me dibuje mediante puntos las distintas aproximaciones que va calculando el método de la bisección, junto con el dibujo de la función; de esa forma conseguimos ver cómo funciona la bisección de una forma visual, pero ya digo, no sé cómo llamar a vectorx y vectory en dibujar
Muchas gracias

[David]

Puedes devolver una lista con los valores...

[vermells]

hola tofol de nuevo porque no pruebas de declarar vectorx y vectory como globales en la funcion biseccion ,es que creo que si las declaras globales puedes usarlas en otras funciones auqnue no estoy muy seguro otra manera seria meter la funcion dibujar en biseccion y otra cosa lo de los indices en vez de i=1 porque no i=0 y te ahorras el vectorx[i-1]

[kawalero]

Aver, para usar el vectorx o vectory dentro de la función dibujar tienes que invocar a biseccion dentro de dibujar y que esta te devuelva esos valores. Es decir,usando un return, igual que en c++.
Lo más comodo y que yo uso mucho para poder poner los datos es hacerme una clase. Ejemplo, una clase que guardará listas:

class datos:
x=[]
y=[]

Y ahora dentro de biseccion haces:

datosbiseccion=datos()

datosbiseccion.x=vectorx
datosbiseccion.y=vectory

return datosbiseccion


Tambien podrias devolverlo directamente:

datosbiseccion=[vectorx,vectory]
return datosbiseccion

Lo de usar una clase es que te servira para cualquier otra función. Si tienes muchas cosas que devolver resulta más cómodo y además puedes ir añadiendo más tipos de datos a la clase conforme avances.

[venkman]

Se me ocurre que una forma un poco más elegante de hacer el método de bisección es hacerlo recursivamente, en lugar de utilizar un bucle.

Ojo, no digo que sea mejor ni peor hacerlo de un modo o de otro. Sólo que a mi me parece un poco más elegante... y vaya, que me apetecía hacerlo xD Así que... he sacado el interprete de Python, le he dado una vueltecilla y he llegado a esto:

Código:

import math

def bisec(f,a,b,fa,fb,delta):
	c = (a+b)/2
	print 'c: ' + repr(c) + ' a: ' + repr(a) + ' b: ' + repr(b) + '\n'
	if (c==a) or (c==b) or math.fabs(a-b)<delta:
		return c,b-a
	fc = f(c)
	if fa*fc<0:
		return bisec(f,a,c,fa,fc,delta)
	else:
		return bisec(f,c,b,fc,fb,delta)

def resolver(f,a,b,delta):
	raiz,precision = bisec(f,a,b,f(a),f(b),delta)
	print 'Raíz encontrada en ' + repr(raiz) + ':\n f('+repr(raiz)+') = ' + repr(f(raiz)) + '\n(precision: '+repr(precision)+')'

Se podría usar con:

Código:

def f(x):
	return x*x*x-x-1

resolver(f,1.0,2.0,1.0e-6)

Por supuesto no he hecho nada de pedir la función al usuario ni nada de eso...