Cómo hacer el método de bisección en Python

Quiero hacer un progtwig Python que ejecute un método de bisección para determinar la raíz de:

f(x) = -26 + 85x - 91x2 +44x3 -8x4 + x5 

El método de Bisección es un método numérico para estimar las raíces de un polinomio f (x).

¿Hay algún pseudocódigo, algoritmos o bibliotecas disponibles que pueda usar para decirme la respuesta?

Aquí hay un código que muestra la técnica básica:

 >>> def samesign(a, b): return a * b > 0 >>> def bisect(func, low, high): 'Find root of continuous function where f(low) and f(high) have opposite signs' assert not samesign(func(low), func(high)) for i in range(54): midpoint = (low + high) / 2.0 if samesign(func(low), func(midpoint)): low = midpoint else: high = midpoint return midpoint >>> def f(x): return -26 + 85*x - 91*x**2 +44*x**3 -8*x**4 + x**5 >>> x = bisect(f, 0, 1) >>> print x, f(x) 0.557025516287 3.74700270811e-16 

Puede ver la solución en una pregunta anterior de Desbordamiento de stack aquí que utiliza scipy.optimize.bisect . O, si su propósito es aprender, el pseudocódigo en la entrada de Wikipedia sobre el método de bisección es una buena guía para hacer su propia implementación en Python, como lo sugiere un comentarista en la pregunta anterior.

Mi implementación es más genérica y, sin embargo, más sencilla que las otras soluciones: (y dominio público)

 def solve(func, x = 0.0, step = 1e3, prec = 1e-10): """Find a root of func(x) using the bisection method. The function may be rising or falling, or a boolean expression, as long as the end points have differing signs or boolean values. Examples: solve(lambda x: x**3 > 1000) to calculate the cubic root of 1000. solve(math.sin, x=6, step=1) to solve sin(x)=0 with x=[6,7). """ test = lambda x: func(x) > 0 # Convert into bool function begin, end = test(x), test(x + step) assert begin is not end # func(x) and func(x+step) must be on opposite sides while abs(step) > prec: step *= 0.5 if test(x + step) is not end: x += step return x 

Con tolerancia:

 # there is only one root def fn(x): return x**3 + 5*x - 9 # define bisection method def bisection( eq, segment, app = 0.3 ): a, b = segment['a'], segment['b'] Fa, Fb = eq(a), eq(b) if Fa * Fb > 0: raise Exception('No change of sign - bisection not possible') while( b - a > app ): x = ( a + b ) / 2.0 f = eq(x) if f * Fa > 0: a = x else: b = x return x #test it print bisection(fn,{'a':0,'b':5}, 0.00003) # => 1.32974624634 

En vivo: http://repl.it/k6q