Is there a function that could be used for calculation of the divergence of the vectorial field? (in matlab) I would expect it exists in numpy/scipy but I can not find it using Google.
I need to calculate
div[A * grad(F)], where
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) is a list of 2D
I know I can calculate divergence like this but do not want to reinvent the wheel. (I would also expect something more optimized) Does anyone have suggestions?
Just a hint for everybody reading that:
the functions above do not compute the divergence of a vector field. they sum the derivatives of a scalar field A:
result = dA/dx + dA/dy
in contrast to a vector field (with three dimensional example):
result = sum dAi/dxi = dAx/dx + dAy/dy + dAz/dz
Vote down for all! It is mathematically simply wrong.
Based on Juh_’s answer, but modified for the correct divergence of a vector field formula
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 uses this exact formula (scroll down to Divergence of a Vector Field)
import numpy as np def divergence(field): "return the divergence of a n-D field" return np.sum(np.gradient(field),axis=0)
The answer of @user2818943 is good, but it can be optimized a little:
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
About 7 times faster:
sum implicitely construct a 3d array from the list of gradient fields which are returned by
np.gradient. This is avoided using
Now, in your question what do you mean by
div[A * grad(F)]?
A * grad(F):
Ais a 2d array, and
grad(f)is a list of 2d arrays. So I considered it means to multiply each gradient field by
- about applying divergence to the (scaled by
A) gradient field is unclear. By definition,
div(F) = d(F)/dx + d(F)/dy + .... I guess this is just an error of formulation.
1, multiplying summed elements
Bi by a same factor
A can be factorized:
Sum(A*Bi) = A*Sum(Bi)
Thus, you can get this weighted gradient simply with:
A is instead a list of factor, one for each dimension, then the solution would be:
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)
What Daniel had modified is the right answer, let me explain self defined func divergence further in more detail :
Function np.gradient() defined as : np.gradient( f) = df/dx, df/dy, df/dz +…
but we need define func divergence as : divergence ( f) = dfx/dx + dfy/dy + dfz/dz +… = np.gradient( fx) + np.gradient(fy) + np.gradient(fz) + …
Let’s test, compare with 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()
Even if it doesn’t have the divergence hand-packaged for you, divergence is pretty simple and the derivative tools they give you in scipy (the ones linked above) give you about 90% of the code prepackaged in a nice, efficient way.
The divergence as a built-in function is included in matlab, but not numpy. This is the sort of thing that it may perhaps be worthwhile to contribute to pylab, an effort to create a viable open-source alternative to matlab.
Edit: Now called http://www.scipy.org/stackspec.html
As far as I can tell, the answer is that there is no native divergence function in numpy. Therefore, the best method for calculating divergence is to sum the components of the gradient vector i.e. calculate the divergence.