【问题标题】：Compute divergence of vector field using python使用python计算向量场的散度
【发布时间】：2012-07-11 06:10:34
【问题描述】：

有没有一个函数可以用来计算矢量场的散度？（matlab）我希望它存在于 numpy/scipy 但我无法使用 Google 找到它。

我需要计算div[A * grad(F)]，其中

F = np.array([[1,2,3,4],[5,6,7,8]]) # (2D numpy ndarray)

A = np.array([[1,2,3,4],[1,2,3,4]]) # (2D numpy ndarray)

所以grad(F) 是一个二维列表ndarrays

我知道我可以像this 那样计算散度，但不想重新发明轮子。（我也希望得到更优化的东西）有人有建议吗？

【问题讨论】：

您需要什么样的订单准确度？你的数组是等间距的吗？
blog.sun.tc/2010/10/jensenshannon-divergence-in-numpy.html
@mgilson 是的，数组是等间距的。我需要双精度。
@ZagorulkinDmitry，Jensen-Shannon 分歧是完全不同的东西
正确的做法，在这里：stackoverflow.com/questions/67970477/…

标签： python numpy scipy

【解决方案1】：

只是给所有阅读者的提示：

上述函数不计算向量场的散度。他们对标量域 A 的导数求和：

结果 = dA/dx + dA/dy

相对于向量场（以三维为例）：

结果 = 总和 dAi/dxi = dAx/dx + dAy/dy + dAz/dz

投反对票！这在数学上完全是错误的。

干杯！

【讨论】：

页面底部有点奇怪。其他答案在数学上确实不正确
这在数学上可能是正确的，但只是迈向答案的第一步。目前的案文没有回答这个问题。下面有更新的答案，实际上回答了手头的问题。
这个答案可以通过告诉我我用谷歌搜索的问题的答案来改进，以及告诉我这里的其他答案不是我想要的。
这里是正确的做法：stackoverflow.com/questions/67970477/…

【解决方案2】：

基于 Juh_ 的回答，但针对向量场公式的正确发散进行了修改

def divergence(f):
    """
    Computes the divergence of the vector field f, corresponding to dFx/dx + dFy/dy + ...
    :param f: List of ndarrays, where every item of the list is one dimension of the vector field
    :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

Matlab's documentation 使用这个精确的公式（向下滚动到向量场的发散）

【讨论】：

【解决方案3】：

import numpy as np

def divergence(field):
    "return the divergence of a n-D field"
    return np.sum(np.gradient(field),axis=0)

【讨论】：

【解决方案4】：

@user2818943 的回答不错，不过可以稍微优化一下：

def divergence(F):
    """ compute the divergence of n-D scalar field `F` """
    return reduce(np.add,np.gradient(F))

时间：

F = np.random.rand(100,100)
timeit reduce(np.add,np.gradient(F))
# 1000 loops, best of 3: 318 us per loop

timeit np.sum(np.gradient(F),axis=0)
# 100 loops, best of 3: 2.27 ms per loop

大约快 7 倍： sum 从np.gradient 返回的梯度场列表中隐式构造一个 3d 数组。这可以避免使用reduce

现在，在您的问题中，div[A * grad(F)] 是什么意思？

关于A * grad(F)：A 是二维数组，grad(f) 是二维数组的列表。所以我认为这意味着将每个梯度场乘以A。
关于将散度应用于（由A 缩放）梯度场尚不清楚。根据定义，div(F) = d(F)/dx + d(F)/dy + ...。我想这只是一个表述错误。

对于1，将求和的元素Bi乘以相同的因子A可以因式分解：

Sum(A*Bi) = A*Sum(Bi)

因此，您可以简单地通过以下方式获得这个加权梯度：A*divergence(F)

如果 ̀A 是一个因子列表，每个维度一个，那么解决方案是：

def weighted_divergence(W,F):
    """
    Return the divergence of n-D array `F` with gradient weighted by `W`

    ̀`W` is a list of factors for each dimension of F: the gradient of `F` over
    the `i`th dimension is multiplied by `W[i]`. Each `W[i]` can be a scalar
    or an array with same (or broadcastable) shape as `F`.
    """
    wGrad = return map(np.multiply, W, np.gradient(F))
    return reduce(np.add,wGrad)

result = weighted_divergence(A,F)

【讨论】：

向量场的散度不是 F = d(Fx)/dx + d(Fy)/dy + ... 吗？正确的公式更像是np.ufunc.reduce(np.add, [np.gradient(F[i], axis=i) for i in range(len(F))])
这一切都取决于F中的数据类型。问题不清楚。我有图像处理经验，因此我认为 F 是 nD 图像，因此梯度是 axis x 然后 y 的导数s（如果有更多）。我总结了一下。如果我理解正确，如果 F 是 m 维的 n*m 2D 向量序列，那么我猜你的公式是正确的。但是，如果 F 大于 2d，我将无法理解

【解决方案5】：

丹尼尔修改的是正确答案，让我更详细地解释自定义函数分歧：

函数np.gradient()定义为：np.gradient(f) = df/dx, df/dy, df/dz +...

但我们需要将 func 散度定义为：散度 (f) = dfx/dx + dfy/dy + dfz/dz +... = np.gradient( fx) + np.gradient(fy) + np.gradient(fz) + ...

我们测试一下，和example of divergence in matlab比较

import numpy as np
import matplotlib.pyplot as plt

NY = 50
ymin = -2.
ymax = 2.
dy = (ymax -ymin )/(NY-1.)

NX = NY
xmin = -2.
xmax = 2.
dx = (xmax -xmin)/(NX-1.)

def divergence(f):
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

y = np.array([ ymin + float(i)*dy for i in range(NY)])
x = np.array([ xmin + float(i)*dx for i in range(NX)])

x, y = np.meshgrid( x, y, indexing = 'ij', sparse = False)

Fx  = np.cos(x + 2*y)
Fy  = np.sin(x - 2*y)

F = [Fx, Fy]
g = divergence(F)

plt.pcolormesh(x, y, g)
plt.colorbar()
plt.savefig( 'Div' + str(NY) +'.png', format = 'png')
plt.show()

--------- 更新版本：包括差异步骤----

感谢@henry 的评论，np.gradient 采用默认步长为1，所以结果可能有些不匹配。我们可以提供自己的差异化步骤。

#https://stackoverflow.com/a/47905007/5845212
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

NY = 50
ymin = -2.
ymax = 2.
dy = (ymax -ymin )/(NY-1.)

NX = NY
xmin = -2.
xmax = 2.
dx = (xmax -xmin)/(NX-1.)


def divergence(f,h):
    """
    div(F) = dFx/dx + dFy/dy + ...
    g = np.gradient(Fx,dx, axis=1)+ np.gradient(Fy,dy, axis=0) #2D
    g = np.gradient(Fx,dx, axis=2)+ np.gradient(Fy,dy, axis=1) +np.gradient(Fz,dz,axis=0) #3D
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], h[i], axis=i) for i in range(num_dims)])

y = np.array([ ymin + float(i)*dy for i in range(NY)])
x = np.array([ xmin + float(i)*dx for i in range(NX)])

x, y = np.meshgrid( x, y, indexing = 'ij', sparse = False)

Fx  = np.cos(x + 2*y)
Fy  = np.sin(x - 2*y)

F = [Fx, Fy]
h = [dx, dy]



print('plotting')
rows = 1
cols = 2
#plt.clf()
plt.figure(figsize=(cols*3.5,rows*3.5))
plt.minorticks_on()


#g = np.gradient(Fx,dx, axis=1)+np.gradient(Fy,dy, axis=0) # equivalent to our func
g = divergence(F,h)
ax = plt.subplot(rows,cols,1,aspect='equal',title='div numerical')
#im=plt.pcolormesh(x, y, g)
im = plt.pcolormesh(x, y, g, shading='nearest', cmap=plt.cm.get_cmap('coolwarm'))
plt.quiver(x,y,Fx,Fy)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im, cax = cax,format='%.1f')


g = -np.sin(x+2*y) -2*np.cos(x-2*y)
ax = plt.subplot(rows,cols,2,aspect='equal',title='div analytical')
im=plt.pcolormesh(x, y, g)
im = plt.pcolormesh(x, y, g, shading='nearest', cmap=plt.cm.get_cmap('coolwarm'))
plt.quiver(x,y,Fx,Fy)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im, cax = cax,format='%.1f')


plt.tight_layout()
plt.savefig( 'divergence.png', format = 'png')
plt.show()

【讨论】：

查看颜色条中的最大数字，它是 0.2，但应该是 3！ stackoverflow.com/questions/67970477/…
@henry 你是对的。实际上，numpy.gradient 颜色条范围比分析结果要小一些。我周围的人通常不关心颜色范围。如果有人向我解释原因，我会很高兴。 @_@
这是因为 no.gradient 使用的步长作为默认值 (1) 与您的步长不对应。见这里：stackoverflow.com/questions/67970477/…
或在这里：stackoverflow.com/questions/67971076/…

【解决方案6】：

基于@paul_chen 的回答，并为 Matplotlib 3.3.0 添加了一些内容（需要传递着色参数，我猜默认颜色图已经改变）

import numpy as np
import matplotlib.pyplot as plt

NY = 20; ymin = -2.; ymax = 2.
dy = (ymax -ymin )/(NY-1.)
NX = NY
xmin = -2.; xmax = 2.
dx = (xmax -xmin)/(NX-1.)

def divergence(f):
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

y = np.array([ ymin + float(i)*dy for i in range(NY)])
x = np.array([ xmin + float(i)*dx for i in range(NX)])

x, y = np.meshgrid( x, y, indexing = 'ij', sparse = False)

Fx  = np.cos(x + 2*y)
Fy  = np.sin(x - 2*y)


F = [Fx, Fy]
g = divergence(F)

plt.pcolormesh(x, y, g, shading='nearest', cmap=plt.cm.get_cmap('coolwarm'))
plt.colorbar()
plt.quiver(x,y,Fx,Fy)
plt.savefig( 'Div.png', format = 'png')

【讨论】：

这并不完全正确！ np。梯度假设步长为 1，但这里的步长不同。这就是为什么您会获得大约 0.6 而不是 3 的最大值。请参见此处：stackoverflow.com/questions/67970477/…

【解决方案7】：

散度作为内置函数包含在matlab中，但不包含numpy。这种事情也许值得为 pylab 做出贡献，努力创建一个可行的开源替代 matlab。

http://wiki.scipy.org/PyLab

编辑：现在称为http://www.scipy.org/stackspec.html

【讨论】：

【解决方案8】：

据我所知，答案是 numpy 中没有原生散度函数。因此，计算散度的最佳方法是将梯度向量的分量相加，即计算散度。

【讨论】：

【解决方案9】：

我不认为@Daniel 的回答是正确的，尤其是当输入是有序的[Fx, Fy, Fz, ...]。

一个简单的测试用例

查看 MATLAB 代码：

a = [1 2 3;1 2 3; 1 2 3];
b = [[7 8 9] ;[1 5 8] ;[2 4 7]];
divergence(a,b)

给出结果：

ans =

   -5.0000   -2.0000         0
   -1.5000   -1.0000         0
    2.0000         0         0

和丹尼尔的解决方案：

def divergence(f):
    """
    Daniel's solution
    Computes the divergence of the vector field f, corresponding to dFx/dx + dFy/dy + ...
    :param f: List of ndarrays, where every item of the list is one dimension of the vector field
    :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])


if __name__ == '__main__':
    a = np.array([[1, 2, 3]] * 3)
    b = np.array([[7, 8, 9], [1, 5, 8], [2, 4, 7]])
    div = divergence([a, b])
    print(div)
    pass

给出：

[[1.  1.  1. ]
 [4.  3.5 3. ]
 [2.  2.5 3. ]]

说明

Daniel 解决方案的错误在于，在 Numpy 中，x 轴是最后一个轴而不是第一个轴。在使用np.gradient(x, axis=0)时，Numpy实际上给出了y方向的梯度（当x为二维数组时）。

我的解决方案

我的解决方案基于丹尼尔的回答。

def divergence(f):
    """
    Computes the divergence of the vector field f, corresponding to dFx/dx + dFy/dy + ...
    :param f: List of ndarrays, where every item of the list is one dimension of the vector field
    :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[num_dims - i - 1], axis=i) for i in range(num_dims)])

在我的测试用例中给出的结果与 MATLAB divergence 相同。

【讨论】：

在笛卡尔索引约定中，MATLAB 遵循 meshgrid 和 NumPy 遵循 meshgrid 中的“xy”索引，x 轴只是第二个轴（x 和 y 交换）。因此 Daniel 的解决方案适用于 [Fy, Fx, Fz, ...]，其中所有 Fn 都在笛卡尔索引中。您的“解决方案”颠倒了所有内容的顺序，并且不适用于维度 > 2，因为它适用于 [...，Fz'，Fy'，Fx']，其中每个 Fn' 的轴顺序相反.对于矩阵索引，Daniel 的解决方案按原样工作。

【解决方案10】：

不知何故，以前计算散度的尝试是错误的！让我告诉你：

我们有以下向量场F：

F(x) = cos(x+2y)
F(y) = sin(x-2y)

如果我们计算散度（使用 Mathematica）：

Div[{Cos[x + 2*y], Sin[x - 2*y]}, {x, y}]

我们得到：

-2 Cos[x - 2 y] - Sin[x + 2 y]

在y [-1,2]和x [-2,2]的范围内有最大值：

N[Max[Table[-2 Cos[x - 2 y] - Sin[x + 2 y], {x, -2, 2 }, {y, -2, 2}]]] = 2.938

使用此处给出的散度方程：

def divergence(f):
        num_dims = len(f)
        return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

我们得到的最大值约为0.625

正确的发散函数：Compute divergence with python

【讨论】：