Compute divergence with python

Question

From this answer, the divergence of a numeric vector field can be computed as such:

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

However, I have noticed that the output seems to depend a lot on the grid resolution, so there seems to be something wrong!

If I look at an example:

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 [-2,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, we get the following plot, for max value vs. resolution (NxN: number of values in x and y-direction). None of these are even close to 3.

Here is the code:

import numpy as np
import matplotlib.pyplot as plt

# Boundaries
ymin = -2.; ymax = 2.
xmin = -2.; xmax = 2.
# Number of points (NxN)
N = 20

# Divergence function
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)])

# Create Meshgrid
x = np.linspace(xmin,xmax, N)
y = np.linspace(ymin,ymax, N)
xx, yy = np.meshgrid(x, y)


# Define 2D Vector Field
Fx  = np.cos(xx + 2*yy)
Fy  = np.sin(xx - 2*yy)

F = np.array([Fx, Fy])
# Compute Divergence
g = divergence(F)

print("Max: ", np.max(g.flatten()))

plt.imshow(g)
plt.colorbar()

Edit: To create the plot:

# %%
a = []
for N in range(20,100):
    # Number of points (NxN)
    # = 20
    # Boundaries
    ymin = -2.; ymax = 2.
    xmin = -2.; xmax = 2.
    
    
    # Deivergence function
    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)])
    
  
    
    # Create Meshgrid
    x = np.linspace(xmin,xmax, N)
    y = np.linspace(ymin,ymax, N)
    xx, yy = np.meshgrid(x, y)
    
    
    # Define 2D Vector Field
    Fx  = np.cos(xx + 2*yy)
    Fy  = np.sin(xx - 2*yy)
    
    F = np.array([Fx, Fy])
    # Compute Divergence
    g = divergence(F)
    
    print("Max: ", np.max(g.flatten()))
    a.append(np.max(g.flatten()))
plt.plot(a)

Answer 1

I realized what the issue was with the help of this answer. The default spacing assumed between two consecutive values in numpy.gradient() is 1. It needs to be changed if there is a different grid.

Hence the divergence function needs to be adapted as such:

Divergence function

def divergence(f,sp):
    """ Computes divergence of vector field 
    f: array -> vector field components [Fx,Fy,Fz,...]
    sp: array -> spacing between points in respecitve directions [spx, spy,spz,...]
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i in range(num_dims)])

Example

a = []
for N in range(20,100):
    # Number of points (NxN)
    # = 20
    # Boundaries
    ymin = -2.; ymax = 2.
    xmin = -2.; xmax = 2.
    
    
    # Divergence function
    def divergence(f,sp):
        num_dims = len(f)
        return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i in range(num_dims)])

    
    # Create Meshgrid
    x = np.linspace(xmin,xmax, N)
    y = np.linspace(ymin,ymax, N)
    xx, yy = np.meshgrid(x, y)
    
    
    # Define 2D Vector Field
    Fx  = np.cos(xx + 2*yy)
    Fy  = np.sin(xx - 2*yy)
    
    F = np.array([Fx, Fy])
    # Compute Divergence
    sp_x = np.diff(x)[0]
    sp_y = np.diff(y)[0]
    sp = [sp_x, sp_y]
    g = divergence(F, sp)
    
    print("Max: ", np.max(g.flatten()))
    a.append(np.max(g.flatten()))
plt.plot(a)

We can see that with increasing resolution, the maximum of the divergence really tends to 3.