Hola, durante las fallas Carlos Sessa me pidió si podía añadir en el código que hice en Python sobre métodos aproximativos de búsqueda de raíces de funciones el método de Newton y le dije que encantado, que no había ningún problema porque el código que hice es totalmente libre, sin ningún tipo de licencia restrictiva, para que todo el mundo aprenda de ello o mejore el código si lo desea. También le prometí publicarlo en el blog, así que eso es lo que voy ha hacer ahora, aunque con algo de retraso

Bueno al meollo. El método de Newton parte de la idea de calcular la raíz (o cero) de una función en un intervalo mediante aproximaciones a la raíz mediante el cálculo de la recta tangente en esas aproximaciones. Para ello utilizamos la ecuación de la recta a la curva de la forma punto-pendiente.

Así pues partiendo de la curva (función) derivable en un punto incial cercano a la raíz que buscamos, sabemos que la ecuación punto-pendiente de la curva f tiene la expresión:

Al ser y todas sus aproximaciones las aproximaciones de la raíz que buscamos, tenemos que , ya que las raíces de cualquier curva, son gráficamente hablando, los puntos de la gráfica de f que cortan en el eje de abcisas, por lo que su ordenada es cero. Así pues tenemos:

Si ahora consideramos que , podemos considerar la expresión anterior en términos de sucesiones. Así pues, sea la sucesión  infinita de términos reales que representan a las distintas aproximaciones de la raíz que buscamos. Como f es derivable en ]a,b[ tenemos que f es continua en [a,b], así pues por el Teorema del Límite Vinculante podemos considerar que existen las sucesiones . Así pues, sustituyendo en la expresión anterior obtenemos:

Si despejamos obtenemos:

Que es la expresión propia del algoritmo de Newton-Raphson si hacemos iterar el número natural "n" desde 0 hasta un número determinado de iteraciones, que vendrá determinado principalmente por la precisión que queramos que tenga la aproximación de nuestra raíz, que es lo que se viene a llamar por el palabro de tolerancia.

Antes es bueno aclarar que una demostración rigurosa de éste método debería hacerse mediante el desarrollo del Polinomio de Taylor de la curva f centrado en el y hasta orden 2, para así obtener una cota del error, pero la omito porque me parece más fácil de entender mediante la ecuación punto-pendiente, que es la que conocen muchos estudiantes de 2º de bachiller.

Hay que recalcar que el método de Newton falla cuando la curva f no es derivable en un entorno pequeño del valor inicial, también falla si en algún punto del entorno del valor inicial su derivada es nula, o próxima a cero hablando desde el punto de vista computacional.

Si queremos buscar entornos del punto inicial para buscar raíces sin dibujar la gráfica de la curva se pueden usar los teoremas de Bolzano y de Rolle. El primero nos da la existencia de la raíz en un intervalo o entorno, y el segundo nos da su unicidad.

Sin más os dejo el código del método de Newton-Raphson escrito en código Python:

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

import Evaluar
 from pylab import *
 from Numeric import *

def newton(df, po, TOL, N):

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

while i<=N:

vectorx[i-1] = po
 Evaluar.dicc_seguro['x'] = po
 f  = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
 df_val = eval(df, {"__builtins__":None}, Evaluar.dicc_seguro)

vectory[i-1] = f

if df_val == 0.0:
 print "La evaluacion de la derivada de la funcion dio 0"
 break

p1 = po - (f/df_val)

if fabs(po-p1) < TOL:
 print "La raiz buscada es: ",po, "con", i-1,"iteraciones"
 break

i += 1
 po = p1

return [vectorx, vectory]

def dibujar(df, po, TOL, N):

x = arange(po-2,po+2,0.1)
 vectores = newton(df, po, TOL, N)

subplot(211)
 plot(x, eval(Evaluar.funcion), linewidth=1.0)
 xlabel('Abcisa')
 ylabel('Ordenada')
 title('Metodo Newton 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()

Os pongo un pantallazo del método usando la función

Bueno os dejo un zip comprimido con la totalidad del código, con todos los métodos que hay hechos hasta ahora.

¡Muchas gracias a Carlos Sessa por su colaboración!
Saludos