Método Newton-Raphson Con Python

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) f:[a,b]\longrightarrow\mathbb{R} derivable en un punto incial x_0 cercano a la raíz que buscamos, sabemos que la ecuación punto-pendiente de la curva f tiene la expresión:

y-y_0 = f\prime (x_0)(x-x_0)

Al ser x_0 y todas sus aproximaciones las aproximaciones de la raíz que buscamos, tenemos que y_0=0 , 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:

y = f\prime (x_0)(x-x_0)

Si ahora consideramos que y=f(x) , podemos considerar la expresión anterior en términos de sucesiones. Así pues, sea x_{n} 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 f(x_n )\quad y\quad f\prime (x_n ) . Así pues, sustituyendo en la expresión anterior obtenemos:

f(x_n ) = f\prime (x_n)(x_{n+1}-x_n)

Si despejamos x_{n+1} obtenemos:

x_{n+1}=x_{n}-\frac{f(x_n)}{f\prime (x_n)}

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

  1. #!/usr/bin/python
  2.  # coding: latin-1
  3.  
  4. import Evaluar
  5.  from pylab import *
  6.  from Numeric import *
  7.  
  8. def newton(df, po, TOL, N):
  9.  
  10. i = 1
  11.  vectorx = zeros(N, Float64)
  12.  vectory = zeros(N, Float64)
  13.  
  14. while i<=N:
  15.  
  16. vectorx[i-1] = po
  17.  Evaluar.dicc_seguro['x'] = po
  18.  f  = eval(Evaluar.funcion, {"__builtins__":None}, Evaluar.dicc_seguro)
  19.  df_val = eval(df, {"__builtins__":None}, Evaluar.dicc_seguro)
  20.  
  21. vectory[i-1] = f
  22.  
  23. if df_val == 0.0:
  24.  print "La evaluacion de la derivada de la funcion dio 0"
  25.  break
  26.  
  27. p1 = po - (f/df_val)
  28.  
  29. if fabs(po-p1) < TOL:
  30.  print "La raiz buscada es: ",po, "con", i-1,"iteraciones"
  31.  break
  32.  
  33. i += 1
  34.  po = p1
  35.  
  36. return [vectorx, vectory]
  37.  
  38. def dibujar(df, po, TOL, N):
  39.  
  40. x = arange(po-2,po+2,0.1)
  41.  vectores = newton(df, po, TOL, N)
  42.  
  43. subplot(211)
  44.  plot(x, eval(Evaluar.funcion), linewidth=1.0)
  45.  xlabel('Abcisa')
  46.  ylabel('Ordenada')
  47.  title('Metodo Newton con f(x)=' + Evaluar.funcion)
  48.  grid(True)
  49.  axhline(linewidth=1, color='r')
  50.  axvline(linewidth=1, color='r')
  51.  
  52. subplot(212)
  53.  plot(vectores[0], vectores[1], 'k.')
  54.  xlabel('Abcisa')
  55.  ylabel('Ordenada')
  56.  grid(True)
  57.  axhline(linewidth=1, color='r')
  58.  axvline(linewidth=1, color='r')
  59.  
  60. show()

Os pongo un pantallazo del método usando la función f(x)=x^2-\pi
Newton Raphson

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

DESCARGAR CODIGO

Yo el código lo ejecuto y edito con Geany (lo tenéis en los repositorios de cualquier distro linuxera importante.
En la Wikipedia tenéis más información sobre éste método->Enlace a la Wikipedia
Cuando tenga algo de tiempo haré el método de Müller, para poder aproximar raíces en el campo de los números complejos.

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

7 Responses to“Método Newton-Raphson Con Python”

  1. 8 abril, 2009 at 0:06 #

    Gracias por el aporte, es muy valioso para mi, estoy empezando a hacer mis pininos con python y en su codigo muestras muy buenos ejemplos de elementos del lenguaje.
    Gracias de nuevo. PD. Visita mi blog: OsoSentado.wordpress.com

  2. 8 abril, 2009 at 2:34 #

    Genial, lo único que no me gusta es el códificado, debería ser utf-8

    Es decir, cambiar:
    # coding: latin-1

    Por:
    # -*- coding: utf-8 -*-

    Saludos

  3. Macarse
    8 abril, 2009 at 7:10 #

    Excelente!

  4. Max
    29 noviembre, 2009 at 18:19 #

    Numeric no es acaso una libreria en desuso ?

    Gracias por tu trabajo!

  5. 13 octubre, 2010 at 6:38 #

    disculpa por que cuando lo ejecuto en geany me apaarece python no se reconoce como lenguaje interno o externo, programa o archivo por lotes ejecutable

  6. Rob
    17 febrero, 2013 at 9:02 #

    Hola tengo una duda como seria la sintaxis correcta para la función: 2x²-x-5=0

  7. 1 marzo, 2013 at 11:13 #

    Sería 2*x**2-x-5

Comments are closed.

Proudly powered by WordPress   Premium Style Theme by www.gopiplus.com
A %d blogueros les gusta esto: