绘制向量函数的最 Pythonic 方式答案

【问题标题】：Most Pythonic way to plot a vector function绘制向量函数的最 Pythonic 方式
【发布时间】：2015-01-20 06:01:07
【问题描述】：

我有一个函数calcField，当给定一个包含两个元素表示一个位置的 numpy 数组时，它返回一个表示该位置电场的数组。要求 matplotlib 为该函数绘制矢量场的最 Pythonic 方式是什么？目前我有这段代码工作，但感觉与 numpy 的精神背道而驰，并且相对不可读。

Y, X = np.mgrid[-3:3:100j, -3:3:100j]

vectors = np.array([[field.calcField(r) for r in row] 
                  for row in [zip(a, b) for a, b in zip(X, Y)]])
U = np.array([[vector[0] for vector in row] for row in vectors])
V = np.array([[vector[1] for vector in row] for row in vectors])

plt.streamplot(X, Y, U, V, color=U, linewidth=2, cmap=plt.cm.autumn)

编辑：根据要求，calcField的代码：

import constants
import numpy as np
import numpy.linalg as l
class Field:
    def __init__(self, charges = []):
       self.charges = charges
    def addCharge(self, charge):
        self.charges = self.charges + [charge]
    def calcField(self, point):
        point = np.array(point)
        return sum([charge.calcField(point) for charge in self.charges])

class PointCharge:
    def __init__(self, q, position):
        self.q = q
        self.position = np.array(position)
    def calcField(self, point):
        return constants.k * self.q * (point - self.position) / l.norm (point - self.position)**3

【问题讨论】：

尝试矢量化函数calcField。一旦掌握了这一点，您就可以摆脱这段代码摘录中的所有其他 for 循环。
补充 Oliver 所说的内容，您可以写 U,V = calcField(X,Y)。只有在值得付出努力的情况下才这样做，即这段代码很慢，或者你知道你将来会经常使用calcField()。
我已经养成了将 N 空间中的向量表示为长度为 N 的 numpy 数组的习惯，因此我可以只处理向量运算，而不是跟踪单个组件。例如，点电荷字段的表达式就是k*q*(pos1 - pos2)/np.linalg.norm(pos1 - pos2)**3。有没有办法在一次对多个输入进行操作的意义上“矢量化”时保持数学的这种可读性？
我们可以看看calcField吗？

标签： python numpy matplotlib

【解决方案1】：

使用流线绘制一组点电荷的电场的代码的矢量化形式可能如下所示：

num_charges = 4
charges = np.random.random_integers(-5,5,num_charges)
charges[charges==0] = 5
charges_positions = np.random.random((num_charges, 2))

y,x = np.mgrid[0:1:40j, 0:1:40j]
xdist = x - charges_positions[:,0].reshape(-1,1,1)
ydist = y - charges_positions[:,1].reshape(-1,1,1)

denom = ((xdist**2 + ydist**2)**1.5)
# Ignoring Coulomb's constant here...
Ex = (charges.reshape(-1,1,1) * xdist / denom).sum(axis=0)
Ey = (charges.reshape(-1,1,1) * ydist / denom).sum(axis=0)

我觉得这比这个替代方案更容易理解，你可能会发现它更具可读性（这是你的问题）：

num_charges = 4
charges = np.random.random_integers(-5,5,(num_charges,1,1))
charges[charges==0] = 5  # only for clarity
positions = np.random.random((2, num_charges,1,1))

y,x = np.mgrid[0:1:40j, 0:1:40j]
M,N = y.shape
xy = np.array([x,y]).reshape(2,1, M,N)
rad_dist = xy - positions
denom = np.linalg.norm(rad_dist, axis=(0))**3

elec_fields = charges * rad_dist / denom
Ex, Ey = elec_fields.sum(axis=1)

你可以很容易地绘制。我将继续使用最后一个代码块的格式（如果您使用第一种格式，则需要交换一些索引）：

pos_charges = charges > 0
neg_charges = charges < 0
f,(ax,ax1) = plt.subplots(1,2)
ax.plot(positions[0, pos_charges], positions[1, pos_charges], 'ro ')
ax.plot(positions[0, neg_charges], positions[1, neg_charges], 'bo ')
ax.streamplot(x,y, Ex, Ey, color='k')
ax.set_aspect('equal', adjustable='box')
ax.set_title('Electric field')
ax.set_xticks([])
ax.set_yticks([])

但是，此时我不再使用类。有时，更轻松地访问矢量化是值得的。在代码的第二种形式中，基本上每个PointCharge.calcField() 的结果都在elec_fields 的第二个轴上，第一个只是这些字段的 x 和 y 分量。

【讨论】：