用python计算矢量场的散度

Question

有没有一个函数可以用来计算矢量场的散度？（在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)是2D ndarray的列表

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

Answer 1

给大家一个提示：

上面的函数不计算向量场的散度。它们求标量场a的导数之和：

结果=dA/dx+dA/dy

与矢量场（以三维为例）相比：

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

投票给所有人！这在数学上是完全错误的。

干杯！

Answer 2

import numpy as np

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

Answer 3

@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)]是什么意思？

1. 关于A * grad(F)：A是二维数组，而grad(f)是二维数组的列表。所以我认为这意味着把每个梯度场乘以A。
2. 关于将散度应用于（按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)