计算和绘制向量场

Question

我正在尝试为一个给定的对象绘制一个潜在场，使用的公式如下：

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

接下来是我使用的代码：

    # 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()

这段代码用来计算一个潜在场。请问我该如何把这个潜在场画得好看一些？另外，给定一个潜在场，最好的方法是什么来把它转换成一个向量场？（向量场是潜在场的负梯度。）

我很感激任何帮助。

我尝试过使用np.gradient()，但结果并不是我预期的那样：

我期待的结果大概是这样的：

编辑：

在代码中更改了两行后：

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

我得到了一个不正确的图像：看起来左右颠倒了（箭头似乎在正确的位置，但场却不对）：

编辑：

通过改成 y, x = np.mgrid[500:-100:200j, -100:1000:200j] 解决了这个问题。你知道为什么吗？

Answer 1

首先，我们在一个普通的网格上进行评估，这和你的示例代码类似。（顺便提一下，你的代码在评估方程时有个错误，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))

接下来，我们需要计算梯度（这是一种简单的有限差分方法，而不是像上面那样通过公式计算导数）：

# 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]))

现在我们可以制作一个“箭头图”，不过结果可能不会像你预期的那样，因为在网格的每个点上都会显示一个箭头：

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

让我们把箭头做得更大。最简单的方法是只绘制每隔n个点的箭头，让matplotlib自动调整大小。这里我们选择每隔3个点。如果你想要更少但更大的箭头，可以把3换成一个更大的整数。

# 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()

这样好一些，但这些箭头还是不太容易看清。更好的可视化方式可能是用一个图像图，并在上面叠加黑色渐变箭头：

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()

理想情况下，我们希望使用不同的颜色映射或改变箭头的颜色。这个部分就留给你自己去尝试吧。你也可以考虑使用等高线图（ax.contour(x, y, p)）或流线图（ax.streamplot(x, y, dx, dy)）。这里给你快速展示一下这些例子：

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()

...为了让效果更炫一些：

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()