Interpolación Lagrange Con Python 3

Hola
En los dos anteriores artículos se vió cómo aproximar una función no polinómica mediante los polinomios de Taylor y de Bernstein. Siguiendo con el tema veamos cómo podemos encontrar un polinomio que pasen por una cantidad finita de puntos. Uno de los primeros polinomios candidatos es el Polinomio De Lagrange.


Veamos cómo podemos calcularlo y dibujarlo en Python 3, utilizaremos: Numpy, Sympy y Matplotlib.
Supongamos un sencillo ejemplo de diferentes contracciones de un resorte dependiendo de las cargas aplicadas.

Carga (Kp) = (5, 10, 15, 20, 25)
Contracción (mm) = (49, 105, 172, 253, 352)

Lo que debemos hacer es con el uso de Numpy declarar 2 matrices fila, una para Carga y otra para Contraccion. Seguidamente, con Sympy hay que utilizar la función interpolating_poly() para calcular el Polinomio de Lagrange. Finalmente, con Matplotlib, dibujamos los puntos de la tabla y el Pol. de Lagrange. Como ves es muy sencillo, sobretodo si tienes práctica con Matlab, Octave, Scilab, etc. De hecho me parece incluso mucho más sencillo con Python, y desde luego no come recursos de tu ordenata. El código es el siguiente:

  1. __author__ = 'tobal'
  2.  
  3. #/usr/bin/env python3
  4. # -*- coding: utf-8 -*-
  5.  
  6. from sympy import *
  7. import matplotlib.pyplot as plt
  8. from  sympy.polys.specialpolys import *
  9. from matplotlib import rcParams
  10. import numpy as np
  11.  
  12. rcParams['text.latex.unicode'] = True
  13. rcParams['text.usetex'] = True
  14. rcParams['text.latex.preamble'] = '\\usepackage{amsthm}', '\\usepackage{amsmath}', '\\usepackage{amssymb}',
  15. '\\usepackage{amsfonts}', '\\usepackage[T1]{fontenc}', '\\usepackage[utf8]{inputenc}'
  16.  
  17. x = symbols('x')
  18.  
  19. init_printing(use_unicode=True)
  20.  
  21. Carga = np.array([5, 10, 15, 20, 25], dtype=float)
  22. Contraccion = np.array([49, 105, 172, 253, 352], dtype=float)
  23. Lagrange = expand(interpolating_poly(len(Carga), x, X=Carga, Y=Contraccion))
  24. pprint(Lagrange)
  25. Lagrange = horner(Lagrange)
  26. P = lambdify(x, Lagrange, 'numpy')
  27.  
  28.  
  29. fig, ax = plt.subplots()
  30. title = r'$Polinomio \; De \; Lagrange \; Grado \; n={}$'.format(len(Carga))
  31. ax.set_title(title)
  32. t = np.linspace(Carga[0], Carga[-1], 1000)
  33. ax.axis([Carga[0] - 1, Carga[-1] + 1, 0, Contraccion[-1] + 48])
  34.      
  35. ax.set_xlabel(r'$Carga\; (Kp)$', color='blue')
  36. ax.set_ylabel(r'$Contraccion\; (mm)$', color='blue')
  37.      
  38. ax.grid('on')
  39.      
  40. ax.plot(t, P(t), 'b-', lw=1.5)
  41. ax.plot(Carga, Contraccion, 'r*', lw=1.5)
  42.      
  43. # Drawing Our Legend
  44. text1 = r'$Pol.\; Lagrange$'
  45. text2 = r'$Puntos\; Tabla$'
  46.      
  47. plt.legend([text1, text2])
  48.      
  49. plt.show()

Dejo imagen del resultado obtenido:

lagrange

En la imagen se aprecia que el resultado es el esperado, y si se me permite, bastante profesional

En fin, un paso más hacia el conocimiento de la programación matemática con Python 3. O mejor dicho, con software totalmente libre, ése software libre al que los ignorantes de éste país le hacen ascos . Pero no nos desviemos del tema, ahora si queremos conocer un valor interpolante en concreto ( o varios ) es muy sencillo. Por ejemplo, si queremos calcular para una carga de 5.75 Kp cuál es su correspondiente Contracción lo podemos hacer añadiendo las siguientes líneas a nuestro código:

  1. pprint('Introduce una carga en Kp a interpolar: ')
  2. CInt = eval(input())
  3. ConInt = Lagrange.subs(x, CInt)
  4. pprint('Para {0} Kp. tenemos una Contracción de {1} m.m. '.format(CInt, ConInt))

El resultado por la Shell es el siguiente:

usuario@nombrepc:~$
Introduce una carga en Kp a interpolar:
5.75
Para 5.75 Kp. tenemos una Contracción de 56.7886773437503 m.m.

Saludos

Pingbacks/Trackbacks

  1. Bitacoras.com - 15 abril, 2014

    Información Bitacoras.com

    Valora en Bitacoras.com: No hay resumen disponible para esta anotación

Comments are closed.

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