13

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 ndarrays

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?

7

10 Answers 10

24

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!

5
  • 5
    Kind of weird how this is at the bottom of the page. The other answers are really mathematically incorrect
    – Daniel
    May 1, 2017 at 18:24
  • 1
    This may be mathematically correct, but is only a first step towards an answer. The current text does not answer the question. There are updated answers below that actually answers the question at hand.
    – vidstige
    Apr 7, 2019 at 17:30
  • 1
    This answer could be improved by also telling me the answer to the question I googled, as well as telling me the other answers here are not what I want.
    – Dast
    Jan 21, 2020 at 17:17
  • Here is the correct way of doing it: stackoverflow.com/questions/67970477/…
    – user7086216
    Jun 14, 2021 at 13:56
  • What??? The Matlab function linked above does indeed calculate the divergence of a vector field. What is that you're referring to, here? It makes no sense whatsoever.
    – Atcold
    Apr 11 at 15:58
12
import numpy as np

def divergence(field):
    "return the divergence of a n-D field"
    return np.sum(np.gradient(field),axis=0)
1
  • field is a scalar field. This does not compute the divergence of a vector field.
    – Atcold
    Apr 11 at 16:00
12

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)

0
11

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

Timeit:

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 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`.
    """
    wGrad = return map(np.multiply, W, np.gradient(F))
    return reduce(np.add,wGrad)

result = weighted_divergence(A,F)
2
  • 4
    Isn't the divergence of a vector field F = d(Fx)/dx + d(Fy)/dy + ... ? The correct formula would be something more like np.ufunc.reduce(np.add, [np.gradient(F[i], axis=i) for i in range(len(F))])
    – Daniel
    May 1, 2017 at 18:23
  • That all depends on the type of data in F. The question is unclear. I have experience in image processing, thus I consider F to be a nD image, and thus the gradient is the derivatives by axis x then y (and more if there are). Which I sum up. If I understand correctly, if F is a n*m 2D sequence of n vector that are m-dimensional, then I guess your formulation is the correct one. HOwever, I would not understand it if F is more than 2d
    – Juh_
    Sep 18, 2017 at 9:27
6

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

enter image description here

---------- UPDATED VERSION: Include the differential Steps----------------

Thank the comment from @henry, the np.gradient take the default step as 1, so the results may have some mismatch. We can provide our own differential steps.

#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()

numerical compare the analytical results

5
  • Look at the maximum number in your colorbar, it's 0.2, but it should be 3! stackoverflow.com/questions/67970477/…
    – henry
    Jun 14, 2021 at 13:55
  • @henry You're right. Indeed, numpy.gradient colorbar range is a bit smaller than analytical results. People around me usually don't care about the color range. I would be happy if someone explains to me why. @_@ Jun 14, 2021 at 21:09
  • It’s because of the step size that no.gradient uses as a default (1) which does not corespondent to your step size. See here: stackoverflow.com/questions/67970477/…
    – henry
    Jun 15, 2021 at 3:52
  • Or here: stackoverflow.com/questions/67971076/…
    – henry
    Jun 15, 2021 at 3:53
  • There's a <typo> and [missing comma] but this post cannot be edited. np.gradient(f) = df/dx, df/dy, df/dz[,] <+>...
    – Atcold
    Apr 11 at 16:13
2

Based on @paul_chen answer, and with some additions for Matplotlib 3.3.0 (a shading param needs to be passed, and default colormap I guess has changed)

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

enter image description here

1
  • This is not entirely correct! np. gradient assumes a step of 1, but here the step is different. This is why you get maximum values that are around 0.6 and not 3. See here: stackoverflow.com/questions/67970477/…
    – henry
    Jun 14, 2021 at 13:54
1

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.

http://wiki.scipy.org/PyLab

Edit: Now called http://www.scipy.org/stackspec.html

1

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.

0
0

I don't think the answer by @Daniel is correct, especially when the input is in order [Fx, Fy, Fz, ...].

A simple test case

See the MATLAB code:

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

which gives the result:

ans =

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

and Daniel's solution:

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

which gives:

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

Explanation

The mistake of Daniel's solution is, in Numpy, the x axis is the last axis instead of the first axis. When using np.gradient(x, axis=0), Numpy actually gives the gradient of y direction (when x is a 2d array).

My solution

There is my solution based on Daniel's answer.

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

which gives the same result as MATLAB divergence in my test case.

1
  • 1
    In the Cartesian indexing convention, which MATLAB follows for meshgrid and NumPy follows for "xy" indexing in meshgrid, the x axis is simply the second axis (x and y are swappped). Thus Daniel's solution will work for [Fy, Fx, Fz, ...], where all Fn are in Cartesian indexing. Your "solution" reverses the order of everything and will not work as intended for dimension > 2, as it works for [..., Fz', Fy', Fx'], where the axes of each Fn' are in opposite order. For matrix indexing, Daniel's solution works as is. Oct 1, 2020 at 14:04
0

Somehow the previous attempts to compute the divergence are wrong! Let me show you:

We have the following vector field F:

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

If we compute the divergence (using Mathematica):

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

we get:

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

which has a maximum value in the range of y [-1,2] and x [-2,2]:

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

Using the divergence equation given here:

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

we get a maximum value of about 0.625

Correct divergence function: Compute divergence with python

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.