Cálculo y dibujo de campos vectoriales.

Estoy tratando de dibujar un campo potencial para un objeto dado usando la siguiente fórmula:

U=-α_goal*e^(-((x-x_goal )^2/a_goal +(y-y_goal^2)/b_goal ) ) 

usando el siguiente código

  # Set limits and number of points in grid xmax = 10.0 xmin = -xmax NX = 20 ymax = 10.0 ymin = -ymax NY = 20 # Make grid and calculate vector components x = linspace(xmin, xmax, NX) y = linspace(ymin, ymax, NY) X, Y = meshgrid(x, y) x_obstacle = 0 y_obstacle = 0 alpha_obstacle = 1 a_obstacle = 1 b_obstacle = 1 P = -alpha_obstacle * exp(-(X - x_obstacle)**2 / a_obstacle + (Y - y_obstacle)**2 / b_obstacle) Ey,Ex = gradient(P) print Ey print Ex QP = quiver(X, Y, Ex, Ey) show() 

Este código calcula un campo potencial. ¿Cómo puedo trazar este campo potencial muy bien? Además, dado un campo potencial, ¿cuál es la mejor manera de convertirlo en un campo vectorial? (campo vectorial es el gradiente negativo del campo potencial).

Apreciaría cualquier ayuda.

He intentado usar np.gradient () pero el resultado no es el que esperaba:

introduzca la descripción de la imagen aquí

Lo que sí espero, es algo a lo largo de estas líneas: introduzca la descripción de la imagen aquí

EDITAR: Después de cambiar las dos líneas en el código:

 y, x = np.mgrid[500:-100:200j, 1000:-100:200j] p = -1 * np.exp(-((x - 893.6)**2 / 1000 + (y - 417.35)**2 / 1000)) 

Tengo un gráfico incorrecto: parece estar invertido de izquierda a derecha (las flechas parecen estar en el lugar correcto pero no en el campo): introduzca la descripción de la imagen aquí EDITAR: Se y, x = np.mgrid[500:-100:200j, -100:1000:200j] cambiando a y, x = np.mgrid[500:-100:200j, -100:1000:200j] Alguna idea de por qué?

En primer lugar, evaluémoslo en una cuadrícula regular, similar a su código de ejemplo. (En una nota al margen, tiene un error en el código para evaluar su ecuación. Falta un negativo dentro de la exp .):

 import numpy as np import matplotlib.pyplot as plt # Set limits and number of points in grid y, x = np.mgrid[10:-10:100j, 10:-10:100j] x_obstacle, y_obstacle = 0.0, 0.0 alpha_obstacle, a_obstacle, b_obstacle = 1.0, 1e3, 2e3 p = -alpha_obstacle * np.exp(-((x - x_obstacle)**2 / a_obstacle + (y - y_obstacle)**2 / b_obstacle)) 

A continuación, necesitaremos calcular el gradiente (esta es una simple diferencia finita, en lugar de calcular analíticamente la derivada de la función anterior):

 # For the absolute values of "dx" and "dy" to mean anything, we'll need to # specify the "cellsize" of our grid. For purely visual purposes, though, # we could get away with just "dy, dx = np.gradient(p)". dy, dx = np.gradient(p, np.diff(y[:2, 0]), np.diff(x[0, :2])) 

Ahora podemos hacer un gráfico de “carcaj”, sin embargo, los resultados probablemente no serán exactamente lo que esperaría, ya que se muestra una flecha en cada punto de la cuadrícula:

 fig, ax = plt.subplots() ax.quiver(x, y, dx, dy, p) ax.set(aspect=1, title='Quiver Plot') plt.show() 

introduzca la descripción de la imagen aquí

Hagamos las flechas más grandes. La forma más sencilla de hacer esto es trazar cada n-ésima flecha y dejar que matplotlib se encargue de la autoescala. Aquí usaremos cada tercer punto. Si desea menos flechas más grandes, cambie el 3 a un número entero más grande.

 # Every 3rd point in each direction. skip = (slice(None, None, 3), slice(None, None, 3)) fig, ax = plt.subplots() ax.quiver(x[skip], y[skip], dx[skip], dy[skip], p[skip]) ax.set(aspect=1, title='Quiver Plot') plt.show() 

introduzca la descripción de la imagen aquí

Mejor, pero esas flechas todavía son bastante difíciles de ver. Una mejor manera de visualizar esto podría ser con un trazado de imagen con flechas negras degradadas superpuestas:

 skip = (slice(None, None, 3), slice(None, None, 3)) fig, ax = plt.subplots() im = ax.imshow(p, extent=[x.min(), x.max(), y.min(), y.max()]) ax.quiver(x[skip], y[skip], dx[skip], dy[skip]) fig.colorbar(im) ax.set(aspect=1, title='Quiver Plot') plt.show() 

introduzca la descripción de la imagen aquí

Lo ideal sería utilizar un mapa de colores diferente o cambiar los colores de las flechas. Te dejo esa parte a ti. También puedes considerar un gráfico de contorno ( ax.contour(x, y, p) ) o un diagtwig de flujo ( ax.streamplot(x, y, dx, dy ). Solo para mostrar un ejemplo rápido de estos:

 fig, ax = plt.subplots() ax.streamplot(x, y, dx, dy, color=p, density=0.5, cmap='gist_earth') cont = ax.contour(x, y, p, cmap='gist_earth') ax.clabel(cont) ax.set(aspect=1, title='Streamplot with contours') plt.show() 

introduzca la descripción de la imagen aquí

… Y por el simple hecho de ser realmente lujoso:

 from matplotlib.patheffects import withStroke fig, ax = plt.subplots() ax.streamplot(x, y, dx, dy, linewidth=500*np.hypot(dx, dy), color=p, density=1.2, cmap='gist_earth') cont = ax.contour(x, y, p, cmap='gist_earth', vmin=p.min(), vmax=p.max()) labels = ax.clabel(cont) plt.setp(labels, path_effects=[withStroke(linewidth=8, foreground='w')]) ax.set(aspect=1, title='Streamplot with contours') plt.show() 

introduzca la descripción de la imagen aquí