Función Meshgrid y gráficos de contorno en Python

Visualizar una gráfica en 2D es relativamente sencillo, puede hacerse a punta de imaginación con algo de práctica. Pero con las funciones de tres variables la situación es más compleja. Sin embargo, siempre es posible obtener información útil a partir de una descripción bidimensional de una figura en 3D gracias a los famosos diagramas de contorno.

Imaginemos una montaña sobre la cual pasamos una serie de planos paralelos a su base a diferentes alturas. La intersección dibujará sobre cada plano la silueta que la montaña posee a esa altura. La proyección de esas siluetas sobre la base de la montaña es una diagrama de contorno. 

Formalmente, dada una función de dos variables $f(x,y)=z$, un diagrama de contorno se define como la curva $f(x,y)=h$ donde $h$ es un valor constante. En otras palabras, los diagramas de contorno son curvas donde la función adopta un valor constate. En el ejemplo de la montaña, si una persona camina a lo largo de una curva de nivel, ni asciende ni desciende. 

Aunque los diagramas de contorno son gráficos 2D, incluso en algunos casos será necesario recurrir a los dispositivos graficadores para crearlos. Vamos a dar un ejemplo. Consideremos una canica que se desliza a lo largo de una alambre parabólico (sin fricción) sometida únicamente a la acción de la gravedad. El Hamiltoniano para este sistema es

$$\displaystyle H(x,p) = \frac{p^2}{2(1+x^2)} + \frac{x^2}{2}$$

siendo $p$ el momentum canónico asociado a la variable de posición $x$. Usando python y la librería matplotlib vamos a crear los diagramas de contorno de esta función. Matplotlib alberga la función contour( ), que recibe tres parámetros relacionados a través de la función $H(x,p)$ los cuales explicaremos después. Antes de esto, notemos que nuestro gráfico de contornos estará sobre el plano $x,p$. Para evaluar la función $H(x,p)$ debemos primero crear una matriz de puntos (o grid) sobre dicho plano. Para esto usamos la función meshgrid( ) de la librería numpy. Esta función recibe dos argumentos que son vectores numpy que contienen información sobre la resolución y el tamaño del grid que deseamos construir. Por ejemplo, si queremos construir una matriz de puntos que tenga un área de 20x20, los argumentos de meshgrid( ) podrían definirse así

	x=np.arange(-10,11,1)
	p=np.arange(-10,11,1)

view raw vectors.py hosted with ❤ by GitHub

En este caso estamos dividiendo los ejes $x$ y $p$ de -10 a 10 en pasos de 1 lo cual define la resolución de nuestro grid. Lo que sigue ahora es crear nuestra matriz de puntos, para lo cual hacemos

X,P=np.meshgrid(x,p)

view raw mesh.py hosted with ❤ by GitHub

Finalmente debemos evaluar la función $H(x,p)$ en cada uno de los puntos de la matriz. Obviamente, nuestra función se evaluará en 400 puntos distintos. La función contour( ) se encarga de hacer esta evaluación automáticamente y de graficar el resultado. Los tres parámetros que recibe esta función son los tres arreglos 400-dimensionales X, Y, h  creados con ayuda de meshgrid( ). Como podemos ver, el arreglo h es precisamente la función $H(x,p)$ evaluada en cada punto del grid. A continuación mostramos el código fuente completo y el resultado.

	import numpy as np
	import matplotlib.pyplot as plt
	x=np.arange(-2, 2, 0.01);
	p=np.arange(-2, 2, 0.01);
	X,P=np.meshgrid(x,p)
	h = (P**2)/(2*(1 + X**2)) + X**2/2.0
	cs = plt.contour(X,P,h)
	plt.clabel(cs, inline=1, fontsize=10)
	plt.show()

view raw contour.py hosted with ❤ by GitHub

Notemos que hay un número que indica el valor constante de la función $H(x,p)$ sobre cada curva. Esto se hace con la función clabel( ) de la línea diez.