【发布时间】：2017-06-26 16:12:33
【问题描述】：
如果我在 matplotlib 中绘制一个矢量场，我通常会明确地为每个组件写下公式，以避免出现形状和广播等问题。然而，在稍微复杂一些的公式中，代码会变得一团糟。
有没有什么方便的方法可以像下面我的（不工作的）伪代码那样更数学地输入涉及向量运算的公式？
# Run with ipython3 notebook
%matplotlib inline
from pylab import *
## The following works, but the mathematical formula is a complete mess to red
def B_dipole(m, a, x,y):
return (3*(x - a[0])*(m[0]*(x - a[0]) + m[1]*(y-a[1]))/((x - a[0])**2 + (y-a[1])**2)**(5/2.0) -m[0]/((x - a[0])**2 + (y-a[1])**2)**(3/2.0),3*(y - a[1])*(m[0]*(x - a[0]) + m[1]*(y-a[1]))/((x - a[0])**2 + (y-a[1])**2)**(5/2.0) -m[1]/((x - a[0])**2 + (y-a[1])**2)**(3/2.0))
## I want something like (but doesn't work)
#def B_dipole(m, a, x,y):
# r = array([x,y])
# rs = r - a ## shifted r
# mrs = dot(m,rs) ## dot product of m and rs
# RS = dot(rs,rs)**(0.5) ## euclidian norm of rs
# ret = 3*mrs*rs/RS**5 - m/RS**3 ## vector/array to return
# return ret
x0, x1=-10,10
y0, y1=-10,10
X=linspace(x0,x1,55)
Y=linspace(y0,y1,55)
X,Y=meshgrid(X, Y)
m = [1,2]
a = [3,4]
Bx,By = B_dipole(m,a,X,Y)
fig = figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1)
ax.streamplot(X, Y, Bx, By,color='black',linewidth=1,density=2)
#ax.quiver(X,Y,Bx,By,color='black',minshaft=2)
show()
输出：
编辑： 我的非工作代码的错误消息：
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-2-43b4694cc590> in <module>()
26 a = [3,4]
27
---> 28 Bx,By = B_dipole(m,a,X,Y)
29
30 fig = figure(figsize=(10,10))
<ipython-input-2-43b4694cc590> in B_dipole(m, a, x, y)
10 def B_dipole(m, a, x,y):
11 r = array([x,y])
---> 12 rs = r - a ## shifted r
13 mrs = dot(m,rs) ## dot product of m and rs
14 RS = dot(rs,rs)**0.5 ## euclidian norm of rs
ValueError: operands could not be broadcast together with shapes (2,55,55) (2,)
如果我不换 r 则会出现错误消息：
--
ValueError Traceback (most recent call last)
<ipython-input-4-e0a352fa4178> in <module>()
23 a = [3,4]
24
---> 25 Bx,By = B_dipole(m,a,X,Y)
26
27 fig = figure(figsize=(10,10))
<ipython-input-4-e0a352fa4178> in B_dipole(m, a, x, y)
8 r = array([x,y])
9 rs = r# - a ## not shifted r
---> 10 mrs = dot(m,rs) ## dot product of m and rs
11 RS = dot(rs,rs)**0.5 ## euclidian norm of rs
12 ret = 3*mrs*rs/RS**5 - m/RS**3 ## vector/array to return
ValueError: shapes (2,) and (2,55,55) not aligned: 2 (dim 0) != 55 (dim 1)
【问题讨论】：
-
您是否尝试使用 numpy 和 numpy.linalg 中的向量运算来精确复制公式？
-
不工作的伪代码正是试图使用这些函数（
dot和abs）。 -
究竟是什么不起作用？有执行错误吗？如果不是，数字将是错误的，因为
RS = abs(rs)计算的是逐项绝对值，而不是欧几里得范数。为此你想要RS = dot(rs,rs)**0.5。 -
为什么
x和y矩阵的大小是55×55？与减去所有条目发生的标量常数相反，当从 55×55 二维向量数组中减去二维向量时会产生歧义，这在所有情况下都无法解析。假设你有一个 2×2 二维向量数组，算法应该如何决定如何分配 3 维的角色？ -
@LutzL：我看到了问题，但我不知道如何规避它。我无法想象任何 matplotlib 的每个熟练用户都可以通过写出所有坐标来输入向量场的公式，就像我在示例中所做的那样。必须有一种方法可以使它更具可读性和可写性，更不容易出错，并且从数学的角度来看更多。
标签： python numpy math matplotlib syntax