Mientras no saco tiempo para algo más elaborado, cálculos de raíces e interpolación con Python! * Cálculo de raíces con el método_de_Newton-Raphson: *

def newton( indice, funcion_der, funcion, x ):

for i in xrange(indice):

 x -= funcion( x ) / float(funcion_der( x ))

return x

El primer parámetro es el número de iteraciones a realizar, el segundo indica la función derivada, después la función, y por último el punto de partida. Por ejemplo: it = 5

print "Raíces [ %s iteraciones]" % it

print "Función: Raíz estimada Valor real"

Newton-Raphson para x

n = newton( it, lambda x: 1, lambda x: x, 10 )

print "x : %s %s" % (n, n)

Newton-Raphson para x + 1

n = newton( it, lambda x: 1, lambda x: x+1, 10 )

print "x + 1 : %s %s" % (n, n+1)

Newton-Raphson para 2x

n = newton( it, lambda x: 2, lambda x: 2*x, 10 )

print "2x : %s %s" % (n, 2n)

Newton-Raphson para x**2

n = newton( it, lambda x: 2x, lambda x: x*2, 10 )

print "x2 : %s %s" % (n, n2)

Salida: Raíces [ 5 iteraciones]

Función: Raíz estimada Valor real

x : 0.0 0.0

x + 1 : -1.0 0.0

2*x : 0.0 0.0

x**2 : 0.3125 0.09765625

* Interpolación con el polinomio_de_Lagrange: *

Este se realiza en dos pasos, primero se obtiene el polinomio de Lagrange para un conjunto de puntos def dd(p1, p2):

return (p1[1] - p2[1]) / ( p1[0][0] - p2[0][-1] )

def lagrange(points):

d = [(points[0][1],[])]

refs = map( lambda x: x[0][0], points)

x = 1

while len(points) > 1:

   p = []

   for i in xrange(len(points) - 1 ):

       p.append((points[i][0] + [points[i+1][0][-1]] ,dd(points[i], points

[i+1])))

   if p[0][1] != 0: # Si se multiplica por 0 ni se considera

       d.append((p[0][1], refs[:x]))



   points = p

   x += 1

return d

Se pasa a lagrange() un array de duplas ( punto, valor en ese punto ).

A continuación se aplica el polinomio a un punto: def xprod( l, x ):

v = 1

for i in l:

   v *= ( x - i )

return v

def solve( l, x ):

return sum(map(lambda y: y[0] * xprod( y[1], x), l)) Se pasa a solve() el polinomio y el punto a esclarecer.

Otra opción es mostrar cómodamente la función:

def fprod( l ):

if len(l) == 0: return ""

a = []

for i in l:

   a.append("x +%s" % -i if i < 0 else "x %s" % -i if i != 0 else "x")

return "(%s)" % ')('.join(a)

def form( l ):

return ' + '.join(map( lambda y: str(y[0]) + fprod( y[1] ),l)) Se le pasaría a form() el polinomio como parámetro.

Por ejemplo: l = lagrange([ # Point, value

     ([1], 1 ),

     ([2], 4 ),

     ([3], 9 )

          ])

print "f(x) = %s"%form(l)

print "\n# Valores de muestra #\n"

for i in xrange( 4, 10):

print "f(%s) = %s" % (i, solve(l, i))

l = lagrange([( [0], 5 )])

print "\n######################\n"

print "f(x) = %s"%form(l)

Salida: f(x) = 1 + 3(x -1) + 1(x -1)*(x -2)

Valores de muestra

f(4) = 16

f(5) = 25

f(6) = 36

f(7) = 49

f(8) = 64

f(9) = 81

f(x) = 5

Lo malo que está quizá demasiado compactado :/

Saludos