计算和绘制矢量场 [英] Computing and drawing vector fields
问题描述
我正在尝试使用以下公式为给定对象绘制势场:
I am trying to draw a potential field for a given object using the following formula:
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()
此代码计算电位场.我怎样才能很好地绘制这个势场?另外,给定一个势场,将其转换为矢量场的最佳方法是什么? (矢量场是势场的负梯度.)
This code calculates a potential field. How can I plot this potential field nicely? Also, given a potential field, what is the best way to convert it to a vector field? (vector field is the minus gradient of the potential field. )
我将不胜感激.
我尝试使用np.gradient(),但结果却不是我所期望的:
I have tried using np.gradient() but the result is not what I have expected:
我期望的是以下几方面的东西:
What I do expect, is something along these lines:
更改代码中的两行之后:
After changing the two lines in the code:
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]
可以解决问题吗?
I have an incorrect plot: it seems to be inverted left-right (arrows seem to be in correct spot but not the field):
Fixed by changing to y, x = np.mgrid[500:-100:200j, -100:1000:200j]
Any idea why?
推荐答案
首先,让我们在常规网格上对其进行评估,类似于您的示例代码. (附带说明,在评估方程式的代码中有错误.在exp
内缺少负数.):
First off, let's evaluate it on a regular grid, similar to your example code. (On a side note, you have an error in the code to evaluate your equation. It's missing a negative inside the 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))
接下来,我们需要计算梯度(与分析计算上述函数的导数相反,这是一个简单的有限差分):
Next, we'll need to calculate the gradient (this is a simple finite-difference, as opposed to analytically calculating the derivative of the function above):
# 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]))
现在我们可以制作一个颤抖"图,但是结果可能不会完全符合您的预期,因为在网格上的每个点都显示了一个箭头:
Now we can make a "quiver" plot, however, the results probably won't be quite what you'd expect, as an arrow is being displayed at every point on the grid:
fig, ax = plt.subplots()
ax.quiver(x, y, dx, dy, p)
ax.set(aspect=1, title='Quiver Plot')
plt.show()
让箭头变大.最简单的方法是绘制第n个箭头,并让matplotlib处理自动缩放.我们将在这里使用每个第3点.如果您想要更少,更大的箭头,请将3更改为更大的整数.
Let's make the arrows bigger. The easiest way to do this is to plot every n-th arrow and let matplotlib handle the autoscaling. We'll use every 3rd point here. If you want fewer, larger arrows, change the 3 to a larger integer number.
# 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()
更好,但是那些箭头仍然很难看到.一种更好的可视化方法可能是覆盖黑色渐变箭头的图像图:
Better, but those arrows are still pretty hard to see. A better way to visualize this might be with an image plot with black gradient arrows overlayed:
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
).只是为了展示一个简单的例子:
Ideally, we'd want to use a different colormap or change the arrow colors. I'll leave that part to you. You might also consider a contour plot (ax.contour(x, y, p)
) or a streamplot (ax.streamplot(x, y, dx, dy
). Just to show a quick example of those:
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()
...而且只是为了变得很花哨:
...And just for the sake of getting really fancy:
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()
这篇关于计算和绘制矢量场的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!