Dibujando Funciones Discontinuas Con Matplotlib

por | 26 enero, 2014

Hola

Estos días he estado investigando sobre dibujar funciones matemáticas discontinuas con Python3 utilizando Matplotlib y Numpy. Me ha costado bastante porque he querido reflejar las gráficas de las funciones en un estilo lo más apropiado para poder utilizarlarlas dando clases.

He tenido 2 problemas principales. El primero ha sido poder situar los ejes coordenados en el centro de la gráfica, en Matplotlib vienen situados en las esquinas por defecto, como en Octave. Googleando he encontrado casi la solución perfecta, sólo faltaba el detalle de que los ejes se extendiesen por toda la imagen aunque no hubiese gráfica dibujada, preguntando en StackOverflow me lo han solucionado en 5 minutos.

El segundo problema ha sido que en las funciones discontinuas la mayoría de programas te dibujan una línea en medio de la discontinuidad, lo cual es erróneo y bastante incómodo. He visto muchas maneras de hacerlo pero sólo me ha funcionado una, que ha sido a la postre la solución más matemática, con lo cual es la mejor. Lo único que me faltaba es conocer el comando correcto a utilizar con Numpy. Enseguida lo explico.

Supongo que muchos se estarán preguntando porqué usar Matplotlib y programar los códigos, cuando herramientas como Mathemathica u otras te lo hacen. Bueno, para empezar Mathemathica será muy bueno pero no es código libre, y muy caro, aunque exista la WolphramAlpha. La otra razón es que en el apartado científico de Python están surgiendo con mucha fuerza alternativas libres como IPython o Sagemath, que vienen a ser poderosos sustitutos de Octave, Scilab, Matlab, Mathemathica, etc… Con IPython y Sagemath puedes trabajar en ciencias tanto con cálculo simbólico como numérico, hay gran variedad de librerías específicas para estadística y con gran cantidad de datos. Con IPython se puede trabajar con clústeres. En ambos programas podemos trabajar desde IDES, terminales enriquecidas o desde el mismo navegador web. Si te pasas por YouTube puedes ver multitud de video-conferencias sobre IPython, Numpy, etc, ninguna española, lógico; aquí siempre vamos los últimos en este tema. Basta de cháchara, empecemos.

Primero veamos una función sencilla definida a trozos continua, para saber cómo dibujar funciones a trozos, la función de Heaviside. Para declarar funciones a trozos utilizaremos la función piecewise de Numpy. El código es sencillo:

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7.  
  8.  
  9. Hv = lambda x: np.piecewise(x, [x < 0.0, x >= 0.0], [0.0, 1.0])
  10.  
  11.  
  12. x = np.linspace(-10, 10, 1000)
  13. plt.axis([x[0], x[-1], -0.1, 1.5])
  14. plt.plot(x, Hv(x), 'b-')
  15. plt.xlabel('x'), plt.ylabel('Hv(x)')
  16. plt.legend(['Hv(x)'])
  17. plt.title('Heaviside Function')
  18. plt.show()

Aquí lo importante es que para que la línea vertical no salga torcida pongamos una buena cantidad de puntos en la malla equiespaciada de x, con la función np.linspace, yo le he puesto 1000. Con eso obtenemos la gráfica siguiente:

Pasemos a una típica función discontinua en un punto y definida a trozos, un trozo es una parábola y el otro una recta.
$latex
f(x) =
\begin{cases}
x^2 & \mbox{si } x< 2 \\ 4 & \mbox{si } x > 2
\end{cases}
$

Aquí el problema es que aunque declaremos la función a trozos con piecewise nos la dibuja totalmente continua, cuando sabemos que es discontinua, la solución fácil es explicitar el punto (2,4) en el comando plot de Matplotlib, que lo dibujamos con la siguiente línea:

  1. ax.plot(x, f(x), 'b-', 2, 4, 'wo', markeredgecolor='b',
  2.         markerfacecolor='w', lw=2.0)

Con markeredgecolor y markerfacecolor conseguimos dibujar el punto como “hueco”, simulando así la discontinuidad. Si queremos dibujar los ejes coordenados en el centro de la gráfica como solemos hacer los docentes debemos hacer el dibujo como un subplot, los comandos para ello son los siguientes:

  1. fig, ax = plt.subplots()
  2. ax.spines['left'].set_position('center')
  3. ax.spines['right'].set_color('none')
  4. ax.spines['bottom'].set_position('center')
  5. ax.spines['top'].set_color('none')
  6. ax.spines['left']
  7. ax.spines['bottom']
  8. ax.xaxis.set_ticks_position('bottom')
  9. ax.yaxis.set_ticks_position('left')

Además, si queremos que los puntos de los ejes salgan bien, sin el problema de que el 0 se repite dos veces, podemos hacerlo con los siguientes comandos:

  1. ticks = []
  2. for i in range(int(x[0]), int(x[-1] + 1), 1):
  3.     ticks.append(i)
  4. ticks.remove(0)
  5. ax.set_xticks(ticks)
  6. ax.set_yticks(ticks)

Lo que hacemos es crearnos una lista de enteros consecutivos llamada ticks, luego borramos el 0 con ticks.remove(0) y con los dos últimos comandos establecemos que nuestros números en ambos ejes van a ser los de nuestra lista ticks. Obviamente, podemos cambiarlos a nuestro gusto fácilmente (reales, caracteres, etc.). El código completo es el siguiente:

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7.  
  8.  
  9. def f(x):
  10.     return np.piecewise(x, [x < 2.0, x > 2.0], [lambda x: x ** 2.0, lambda x: 4.0])
  11.  
  12.  
  13. fig, ax = plt.subplots()
  14. x = np.linspace(-5.0, 5.0, 1000)
  15. ax.axis([x[0], x[-1], x[0], x[-1]])
  16. ax.spines['left'].set_position('center')
  17. ax.spines['right'].set_color('none')
  18. ax.spines['bottom'].set_position('center')
  19. ax.spines['top'].set_color('none')
  20. ax.spines['left']
  21. ax.spines['bottom']
  22. ax.xaxis.set_ticks_position('bottom')
  23. ax.yaxis.set_ticks_position('left')
  24. ticks = []
  25. for i in range(int(x[0]), int(x[-1] + 1), 1):
  26.     ticks.append(i)
  27. ticks.remove(0)
  28. ax.set_xticks(ticks)
  29. ax.set_yticks(ticks)
  30. ax.plot(x, f(x), 'b-', 2, 4, 'wo', markeredgecolor='b',
  31.         markerfacecolor='w', lw=2.0)
  32. ax.legend([r'$f(x)$'], loc='lower right')
  33. ax.set_title(r'$Funci\'on\; A\; Trozos$')
  34. ax.grid('on')
  35. plt.show()

Pasemos ahora a una función racional, f(x)=\frac{x^3}{x^2-x-6} , al dibujarla con Matplotlib tenemos que nos dibuja las asíntotas verticales como líneas continuas, es decir, Matplotlib tiene la costumbre de convertir las funciones discontinuas con polos en funciones continuas; supongo que se debe a que Matplotlib dibuja uniendo puntos si en plot le decimos que el estilo sea línea continua. Para evitarlo:

  1. x = np.linspace(-15.0, 15.0, 1000)
  2. pos = np.where(np.abs(np.diff(f(x))) >= 10.0)[0]
  3. x = np.insert(x, pos, np.nan)

La solución esta utilizar Numpy, definimos los x como una malla equiespaciada con linspace, seguidamente creamos un posicional, pos, que lo que hace es buscar aquellos x de la malla tales que el módulo de la diferencia entre dos f(x) consecutivos sea mayor que, igual o menor que una condición. Es decir; intentamos medir la distancia que hay entre las imágenes de dos puntos consecutivos de nuestra función, vamos que medimos si el salto de discontinuidad en nuestra gráfica. Seguidamente redefinimos x insertando los pos anteriores como tipo numérico np.nan de Numpy (NAN = Not A Number), de esa forma evitamos que dibuje las discontinuidades. Esto no lo haremos siempre, sólo en los casos en que nos dibuje las discontinuidades, y adaptándolo a nuestra función. El código completo es:

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7.  
  8.  
  9. f = lambda x: pow(x, 3) / (pow(x, 2) - x - 6)
  10.  
  11.  
  12. fig, ax = plt.subplots()
  13. x = np.linspace(-15.0, 15.0, 1000)
  14. pos = np.where(np.abs(np.diff(f(x))) >= 10.0)[0]
  15. x = np.insert(x, pos, np.nan)
  16. ax.axis([x[0], x[-1], -15.5, 15.5])
  17. ax.spines['left'].set_position('center')
  18. ax.spines['right'].set_color('none')
  19. ax.spines['bottom'].set_position('center')
  20. ax.spines['top'].set_color('none')
  21. ax.spines['left']
  22. ax.spines['bottom']
  23. ax.xaxis.set_ticks_position('bottom')
  24. ax.yaxis.set_ticks_position('left')
  25. ticks = []
  26. for i in range(-15, 16, 5):
  27.     ticks.append(i)
  28. ticks.remove(0)
  29. ax.set_xticks(ticks)
  30. ax.set_yticks(ticks)
  31. ax.plot(x, f(x), color='b', linestyle='-', lw=1.5)
  32. ax.plot(x, x + 1.0, color='r', linestyle='--', lw=2.0)
  33. ax.axvline(x=-2.0, ymin=-15.0, ymax=15.0, linewidth=2.0,
  34.            color='g', linestyle='--')
  35. ax.axvline(x=3.0, ymin=-15.0, ymax=15.0, linewidth=2.0,
  36.            color='brown', linestyle='--')
  37. ax.legend([r'$f(x)=\frac{x^3}{x^2-x-6}$', r'$y=x+1$', r'$x=-2$', r'$x=3$'], loc='lower right')
  38. ax.annotate(r'$OX$', xy=(13.5, 0.75), size=16, color='black')
  39. ax.annotate(r'$OY$', xy=(0.25, 14.0), size=16, color='black')
  40. ax.annotate(r'$x = 3$', xy=(3.25, 1), size=16, color='brown')
  41. ax.annotate(r'$x = -2$', xy=(-4.5, 1), size= 16, color='g')
  42. ax.annotate(r'$y = x + 1$', xy=(-10, -5), rotation=35, size=16, color='r')
  43. ax.set_title(r'$Funci\'on\; Discontinua$', fontsize=18)
  44. ax.grid('on')
  45. plt.show()

Pasemos a otra función típica, la función piso, en esta también tendremos que hacer lo de borrar las discontinuidades, como antes. Aquí sabemos que los puntos extremos de cada segmento de la gráfica son continuidades por la izqda. y discontinuidades por la derecha. Para reflejar esto de la forma típica que conocemos lo haremos con los siguientes comandos:

  1. igual = np.arange(-6, 6, 1)
  2. igual_mas_uno = igual + np.ones(12, np.int)
  3. ax.plot(igual, igual, 'bo', markeredgecolor='b', markerfacecolor='b',
  4.         lw=2.0,  label='_nolegend_')
  5. ax.plot(igual_mas_uno, igual, 'bo', markeredgecolor='b', markerfacecolor='w',
  6.         lw=2.0,  label='_nolegend_')

Lo que hacemos es crearnos dos arrays (igual, igual_mas_uno) para dibujar con plot los puntos. Al dibujarlos, si en el código ponemos una leyenda Matplotlib nos añadirá los puntos como funciones, para evitarlo añadimos en plot la opción label=’_nolegend_’. El código:

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7.  
  8.  
  9. fig, ax = plt.subplots()
  10. x = np.linspace(-6.0, 6.0, 1000)
  11. pos = np.where(np.abs(np.diff(np.floor(x))) >= 1.0)[0] + 1
  12. x = np.insert(x, pos, np.nan)
  13. ax.axis([x[0] - 0.5, x[-1] + 0.5, x[0] - 0.5, x[-1] + 0.5])
  14. ax.spines['left'].set_position('center')
  15. ax.spines['right'].set_color('none')
  16. ax.spines['bottom'].set_position('center')
  17. ax.spines['top'].set_color('none')
  18. ax.spines['left']
  19. ax.spines['bottom']
  20. ax.xaxis.set_ticks_position('bottom')
  21. ax.yaxis.set_ticks_position('left')
  22. ticks = []
  23. for i in range(int(x[0]), int(x[-1] + 1), 1):
  24.     ticks.append(i)
  25. ticks.remove(0)
  26. ax.set_xticks(ticks)
  27. ax.set_yticks(ticks)
  28. ax.plot(x, np.floor(x), color='b', linestyle='-', lw=2.0)
  29. igual = np.arange(-6, 6, 1)
  30. igual_mas_uno = igual + np.ones(12, np.int)
  31. ax.plot(igual, igual, 'bo', markeredgecolor='b', markerfacecolor='b',
  32.         lw=2.0,  label='_nolegend_')
  33. ax.plot(igual_mas_uno, igual, 'bo', markeredgecolor='b', markerfacecolor='w',
  34.         lw=2.0,  label='_nolegend_')
  35. ax.legend([r'$f(x)=\lfloor x \rfloor$'], loc='lower right')
  36. ax.annotate(r'$OX$', xy=(x[-1] - .5, 0.25), size=16, color='black')
  37. ax.annotate(r'$OY$', xy=(0.25, x[-1]), size=16, color='black')
  38. ax.set_title(r'$Funci\'on\; Piso$', fontsize=18)
  39. ax.grid('on')
  40. plt.show()

De forma similar podemos dibujar la función techo:

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7.  
  8.  
  9. fig, ax = plt.subplots()
  10. x = np.linspace(-6.0, 6.0, 1000)
  11. pos = np.where(np.abs(np.diff(np.ceil(x))) >= 1.0)[0] + 1
  12. x = np.insert(x, pos, np.nan)
  13. ax.axis([x[0] - 0.5, x[-1] + 0.5, x[0] - 0.5, x[-1] + 0.5])
  14. ax.spines['left'].set_position('center')
  15. ax.spines['right'].set_color('none')
  16. ax.spines['bottom'].set_position('center')
  17. ax.spines['top'].set_color('none')
  18. ax.spines['left']
  19. ax.spines['bottom']
  20. ax.xaxis.set_ticks_position('bottom')
  21. ax.yaxis.set_ticks_position('left')
  22. ticks = []
  23. for i in range(int(x[0]), int(x[-1] + 1), 1):
  24.     ticks.append(i)
  25. ticks.remove(0)
  26. ax.set_xticks(ticks)
  27. ax.set_yticks(ticks)
  28. ax.plot(x, np.ceil(x), color='b', linestyle='-', lw=2.0)
  29. igual = np.arange(-5, 7, 1)
  30. igual_mas_uno = igual - np.ones(12, np.int)
  31. ax.plot(igual, igual, 'bo', markeredgecolor='b', markerfacecolor='b',
  32.         lw=2.0,  label='_nolegend_')
  33. ax.plot(igual_mas_uno, igual, 'bo', markeredgecolor='b', markerfacecolor='w',
  34.         lw=2.0,  label='_nolegend_')
  35. ax.legend([r'$f(x)=\lceil x \rceil$'], loc='lower right')
  36. ax.annotate(r'$OX$', xy=(x[-1] - .5, 0.25), size=16, color='black')
  37. ax.annotate(r'$OY$', xy=(0.25, x[-1]), size=16, color='black')
  38. ax.set_title(r'$Funci\'on\; Techo$', fontsize=18)
  39. ax.grid('on')
  40. plt.show()

Finalicemos con el código f(x)=x-\lceil x \rceil :

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6.  
  7.  
  8. fig, ax = plt.subplots()
  9. x = np.linspace(-6.0, 6.0, 1000)
  10. pos = np.where(np.abs(np.diff(np.ceil(x))) == 1.0)[0] + 1
  11. x = np.insert(x, pos, np.nan)
  12. ax.axis([x[0] - 0.5, x[-1] + 0.5, x[0] - 0.5, x[-1] + 0.5])
  13. ax.spines['left'].set_position('center')
  14. ax.spines['right'].set_color('none')
  15. ax.spines['bottom'].set_position('center')
  16. ax.spines['top'].set_color('none')
  17. ax.spines['left']
  18. ax.spines['bottom'].set_smart_bounds(True)
  19. ax.xaxis.set_ticks_position('bottom')
  20. ax.yaxis.set_ticks_position('left')
  21. ticks = []
  22. for i in range(int(x[0]), int(x[-1] + 1), 1):
  23.     ticks.append(i)
  24. ticks.remove(0)
  25. ax.set_xticks(ticks)
  26. ax.set_yticks(ticks)
  27. ax.plot(x, x - np.ceil(x), color='b', linestyle='-', lw=2.0)
  28. igual = np.arange(-5, 7, 1)
  29. igual2 = np.arange(-6, 6, 1)
  30. ax.plot(igual, np.zeros(12, np.int), 'bo', markeredgecolor='b', markerfacecolor='b',
  31.         lw=2.0,  label='_nolegend_')
  32. ax.plot(igual2, -1 * np.ones(12, np.int), 'bo', markeredgecolor='b', markerfacecolor='w',
  33.         lw=2.0,  label='_nolegend_')
  34. ax.legend([r'$f(x)=x-\lceil x \rceil$'], loc='lower right')
  35. ax.annotate(r'$OX$', xy=(x[-1] - .5, 0.25), size=16, color='black')
  36. ax.annotate(r'$OY$', xy=(0.25, x[-1]), size=16, color='black')
  37. ax.set_title(r'$Funci\'on\; f(x)=x-\lceil x \rceil$', fontsize=18)
  38. ax.grid('on')
  39. plt.show()

Con todos estos códigos podemos arreglarnos para poder dibujar infinidad de funciones discontinuas con Matplotlib y Numpy en Python 3, porque obtenemos experiencia en funciones difíciles de dibujar por peculiaridades del software. Además, hemos aprendido a dibujarlas con buen estilo y adaptadas a la educación secundaria, solventando los típicos problemas que nos surgen al querer cambiar la configuración por defecto del software.

Saludos

Un pensamiento en “Dibujando Funciones Discontinuas Con Matplotlib

  1. Pingback: Bitacoras.com

Los comentarios están cerrados.