Compute divergence of vector field using python

Posted on

Question :

Compute divergence of vector field using python

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

so `grad(F)` is a list of 2D `ndarray`s

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.

Cheers!

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"
``````

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` """
``````

Timeit:

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

# 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 `reduce`

Now, in your question what do you mean by `div[A * grad(F)]`?

1. about `A * grad(F)`: `A` is a 2d array, and `grad(f)` is a list of 2d arrays. So I considered it means to multiply each gradient field by `A`.
2. 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.

For `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*divergence(F)`

If ?`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`.
"""

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.