12

ベクトル場の発散の計算に使用できる関数はありますか? ( 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)2Dndarrayのリストも同様です

このように発散を計算できることはわかっていますが、車輪の再発明はしたくありません。（私はもっと最適化されたものも期待しています）誰か提案がありますか？

4

10 に答える 10

10

@ 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によって返される勾配フィールドのリストから暗黙的に 3 次元配列を作成しnp.gradientます。これは使用して回避されますreduce

さて、あなたの質問では、 とはどういう意味 div[A * grad(F)]ですか?

  1. about A * grad(F):Aは 2 次元配列であり、2 次元配列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)

À が各次元に 1 つずつ、因子のリストである場合、解は次のようになります。

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)
于 2014-01-15T10:01:59.573 に答える
0

特に入力が順番にある場合、@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. ]]

説明

ダニエルの解決策の間違いは、Numpy では、x 軸が最初の軸ではなく最後の軸であることです。を使用するnp.gradient(x, axis=0)と、Numpy は実際にy 方向の勾配を与えます(x が 2 次元配列の場合)。

私の解決策

ダニエルの答えに基づく私の解決策があります。 

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)])

divergence私のテスト ケースでは、MATLAB と同じ結果が得られます。

于 2019-05-05T05:49:20.870 に答える